domain: sdu spi.signal stat

text:
Je Remercie

Christophe Baehr

Keywords:

La qualité des informations produites par les Systèmes de Navigation tel que GPS / EGNOS dépendent de beaucoup de facteurs. Nous nous intéressons plus précisément aux évènements ionosphériques qui sont les plus importants dans la déficience de ces services. Cette étude nous place en amont de ces évènements, nous étudions la relation qui peut exister entre l'activité solaire et géomagnétique. Les évènements ionosphériques dépendent directement du géomagnétisme mais nous n'avons pas encore une idée précise de la relation avec le soleil. Par l'utilisation d'outils de traitement de signal ainsi que par optimisation de modèle, nous prouvons que les activités géomagnétique et solaire sont simplement liées par une relation linéaire. Cette relation évolue dans le temps mais reste de la même nature. Nous avons également procédé à une étude des différents cycles que contiennent les données solaires et géomagnétiques et il s'avère que certains soient communs entre les deux activités. Enfin nous faisons l'étude des occurrences des orages géomagnétique que nous considérons comme des évènements rares cependant il est possible de leur attribuer une probabilité selon leur puissance. Mots-clés :

Table des figures

Introduction

Le positionnement par GPS en navigation est extrêmement complexe à mettre en place.

Pour repérer un point sur le globe terrestre, il faut mobiliser quatre satellites : trois pour la tri- Ces systèmes répondent à certains besoins que l'on mesure au travers des Performances GNSS au nombre de quatre, la précision, l'intégrité, la continuité et la disponibilité.

-La précision concerne l'estimation des positions horizontales, verticales mais aussi dans le temps.

-L'intégrité correspond au degré de confiance que nous pouvons accorder aux données de position obtenues.

-La continuité : cette performance permet d'avoir accès au positionnement intègre sans interruption.

-La disponibilité nous permet d'utiliser le service à n'importe quel moment.

Pour assurer l'intégrité des données, chaque pays ou région développé possède un système de correction qui permet d'avoir des positionnements qui se rapprochent le plus de la réalité.

L

Bruit récepteur 1m

Multi-trajets 1m

Disponibilité sélective 0 (24m)

Total UERE 5m

Dilution moyenne de la précision (HDOP) 1,5

Erreur de précision horizontale (HPE 95%) 15m(=5 × 1, 5 × 2)

formations nécessaires à ces calculs de correction peuvent devenir manquantes impactant par la même occasion la disponibilité et continuité du service. Toutes ces erreurs ont des conséquences sur l'efficacité de la navigation, celles liées à l'activité ionosphérique sont les plus importantes.

Il y a donc de réels enjeux à connaître finement les différents évènements géomagnétiques et les relations qu'ils entretiennent avec l'activité solaire. Le tableau 1 nous donne un ordre d'idée des erreurs qui sont provoquées par les différentes causes.

La problématique

Cycles et occurrences

Les hypothèses d'un cycle de 27 jours dans les indices semblent être fondées (Soukharev, 2000), une étude a permi d'examiner les amplitudes du cycle de 27 jours sur les indices Ap, F10.7cm et FoF2 [START_REF] Apostolov | The 22-year cycle in the geomagnetic 27-day recurrences reflecting on the F2-layer ionization[END_REF] Une correction est ensuite appliquée sur ces données pour les erreurs liées à l'absorption atmosphérique. Le F10.7 s'exprime en SFU (Solar Flux Unit).

Indices géomagnétiques

Pour l'activité géomagnétique, nous avons de nombreux indices (figure 1). Toutes ces don-

Mise en forme des indices

f (x) = ∑ n∈Z < f , e n > e n = 1 √ 2π ∑ n∈Z < f , e n > e i.n.x (2)
Où < f , e n > correspond à la contribution de chaque fréquence et aux harmoniques/fréquences du signal.

Le signal étant transformé, nous avons la possibilité de le traiter sous forme de fréquences.

Les coefficients de la transformée de Fourier appartiennent au domaine des exponentielles complexes. La sélection des fréquences est réalisée dans ce sous ensemble, et après ce traitement nous pouvons utiliser l'inverse de la transformée de Fourier (IFFT) pour recomposer le signal traité.

La

R x,x (τ) = E[x(t)x(t + τ)] (3)
pour tout t avec τ la durée entre deux instants de mesure.

S x,x (ω) = +∞ -∞ R x,x (t)e -2iπωt dt (
ψ ab (t) = 1 |a| ψ t -b a (5)
où a > 0 est le facteur d'espace et b le facteur temporel [START_REF] Meyer | Sites internet • Sites de références bibliographiques[END_REF].

(T wav f )(a, b) = R f (t)ψ t -b a dt (6) f (x) = 1 C ψ R R (T wav f )(a, b)ψ ab (x) 1 a 2 da.db (7)
où C représente l'intégrale de la FFT sur l'ondelette mère.

Spline naturelle

Occurrences des évènements rares

Etude des différents cycles

Pour cette partie nous nous attachons à étudier de manière plus fine les différents cycles qui peuvent être présents dans les signaux que nous avons choisis. La méthode utilisée pour cette analyse est principalement la transformée de Fourier, bien adaptée pour l'étude des signaux sur des plages de fréquences variées.

L'analyse est réalisée sur les densités spectrales de puissance des deux signaux. Il est important de les comparer entre elles (figure 4), afin de pouvoir déceler les fréquences communes où nous avons des pics de puissance. Nous pouvons d'ors et déjà constater deux plages de fréquences où les deux indices présentent de forts pics de puissance : entre 1.10 -4 et 3.10 -4 jours -1

(11ans) et entre 3.10 -2 et 4.10 -2 jours -1 (27jours). Ensuite l'indice Aa présente d'autres pics de puissances absents sur la DSP des taches solaires. Les deux autres pics principaux se situent entre 5.10 -3 et 6.10 -3 jours -1 (6mois) et entre 7.10 -2 et 8.10 -2 jours -1 (13jours (Mendoza, 1996 ;[START_REF] Hathaway | A synthesis of solar cycle prediction techniques[END_REF][START_REF] Clúa De Gonzalez | Annual variation of geomagnetic activity[END_REF][START_REF] Ahluwalia | Preliminary forecast for the peak of solar activity cycle 24[END_REF]Prestes et al., 2005), de plus ces deux indices présentent une très grande plage de recensement (entre 1868 et 2011, en commun). L'analyse des corrélations entre les indices géomagnétiques montre des coefficients très proches de 1 pour la majorité. Seul le Dst s'illustre par des coefficients toujours aux environs de 0,70 vis à vis des autres indices géomagnétiques. Cet indice est très particulier car sa mesure est effectuée seulement au niveau de l'équateur contrairement aux autres. De plus il est souvent utilisé dans le but de repérer les orages magnétiques dans les études d'occurrences (Bravo et al., 1997 ;[START_REF] Dubey | Solar wind streams and intense geomagnetic storms observed during 1986-1991[END_REF][START_REF] Echer | Statistical studies of geomagnetic storms with peak Dst ≤ -50nT from 1957 to 2008[END_REF] avec le Ap et le Kp.

• L'analyse spectrale

Nous avons pu mettre en évidence plusieurs cycles réguliers dans les signaux bruts de SSN et Aa. Les deux cycles communs (pour lesquels le géomagnétisme est directement lié à l'activité solaire), sont le cycle d'environ 11 ans correspondant au cycle solaire, et le cycle de 27 jours qui est dû à la rotation du soleil dont la révolution se fait au bout de 27 jours. Un cycle de 13 jours est aussi visible sur l'indice Aa, comme le cycle de 5,5ans pour celui de 11ans, il serait une harmonique du cycle de 27 jours causé par l'activité solaire. Ensuite nous avons pu définir un cycle saisonnier à forte puissance sur le spectre de l'indice géomagnétique choisi n'apparaissant pas sur l'activité solaire, ce qui semble jusqu'ici normal. Nous avons ensuite quelques pics de puissance correspondant à des cycles d'environ 17 mois et 22 mois. D'autres analyses spectrales (Prestes et al., 2005 ;[START_REF] Střeštík | Variations in solar and geomagnetic activity with periods near to 1.3 year[END_REF] ont pu mettre en avant des cycles que nous n'avons pas retrouvés lors de nos travaux.

• Etude de la relation Adler, 1996 ;[START_REF] Hathaway | A synthesis of solar cycle prediction techniques[END_REF][START_REF] Verbanac | Four decades of geomagnetic and solar activity : 1960-2001[END_REF][START_REF] Ahluwalia | Solar cycle 23 prediction update[END_REF]. De plus nous avons pu démontrer un déphasage entre les signaux, la sinusoïde géomagnétique semble être en retard par rapport à la sinusoïde solaire, en recalant les signaux nous avons fait correspondre une activité solaire avec l'activité géomagnétique répondant environ 16mois après. Certains travaux [START_REF] Ahluwalia | Solar cycle 23 prediction update[END_REF][START_REF] Verbanac | Four decades of geomagnetic and solar activity : 1960-2001[END_REF]

• Conclusion

Au cours de cette étude nous avons pu définir la relation qui existe entre l'activité solaire et géomagnétique, le modèle que nous avons réalisé permet d'obtenir une prédiction de la tendance géomagnétique des 16 prochains mois. En premier lieu, l'étude des cycles a permis de mettre en avant les points communs entre les deux signaux à l'aide d'une analyse spectrale. L'étude des occurrences a permis d'obtenir une relation entre la probabilité d'occurrence des orages et l'activité moyenne de chaque cycle solaire. Il est toujours possible de compléter cette dernière étude. Les méthodes pour effectuer la prédiction peuvent être mises en oeuvre, il manque juste les bonnes variables pouvant expliquer l'occurrence des orages. Ces travaux pourront être réutilisés dans le but de prévenir les dysfonctionnements des systèmes GPS dus aux évènements ionosphériques. Une mise à jour de la prédiction de la tendance géomagnétique est possible par l'ajout des dernières données issues des mesures de l'activité solaire. Cette étude est novatrice par rapport aux objectifs que nous visons. Elle permet d'avoir un aperçu du niveau d'activité géomagnétique des prochains mois, de connaître le risque d'avoir un orage magnétique sévère, tout ceci dans le but de préserver les systèmes GPS sans lesquels la navigation serait impossible.

2. 1

 1 Présentation de l'organisme d'accueil . . . . . . . . . . . . . . . . . . . . . . 2.2 Etat de l'art du GNSS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 La problématique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Recherche bibliographique . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Prédictions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.3 Cycles et occurrences . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Objectifs du stage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Matériel & Méthodes 3.1 Les différents indices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Indices solaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Indices géomagnétiques . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Mise en forme des indices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Généralités sur le traitement du signal . . . . . . . . . . . . . . . . . . . . . . 3.4 Les méthodes de traitement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 La Transformée de Fourier Rapide (FFT) . . . . . . . . . . . . . . . . 3.4.2 La Densité Spectrale des Puissances (DSP) . . . . . . . . . . . . . . . 3.4.3 Transformée par ondelettes . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Spline naturelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Optimisation, régression non-linéaire . . . . . . . . . . . . . . . . . . . . . . . 3.7 Occurrences des évènements rares . . . . . . . . . . . . . . . . . . . . . . . . 4 Résultats 4.1 Les relations entre les différents indices géomagnétiques . . . . . . . . . . . . 4.2 Etude des différents cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Mise en évidence d'une relation activité solaire/activité géomagnétique . . . . 4.3.1 Traitement par ondelette de Daubechies . . . . . . . . . . . . . . . . . 4.3.2 Traitement par ondelette de Morlet . . . . . . . . . . . . . . . . . . . . 4.4 Analyse de phase entre les signaux . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Définition de la relation entre SSN et Aa . . . . . . . . . . . . . . . . . . . . . 4.6 Occurrence des maximums . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6.1 Analyse sur l'ensemble des données . . . . . . . . . . . . . . . . . . . 4.6.2 Etude de chaque cycle . . . . . . . . . . . . . . . . . . . . . . . . . .

1234567923

  Ensemble des indices solaires et géomagnétiques en fonction du temps. . . . . Traitement de l'indice Aa par la transformée de Fourier . . . . . . . . . . . . . Signaux recomposés de tous les indices . . . . . . . . . . . . . . . . . . . . . Spectre des Puissances des signaux SSN et Aa . . . . . . . . . . . . . . . . . . Graphiques représentant les indices SSN et Aa mensuels en fonction du temps. . Décomposition multi-échelle des signaux par l'ondelette « db5 » . . . . . . . . Image des transformations multi-échelles des signaux avec l'ondelette de Morlet 8 Graphiques de corrélation entre SSN et Aa sur l'ensemble des échelles . . . . . Graphique des signaux recomposés après ondelette de Daubechies 5 . . . . . . 10 Spectre de déphasage entre l'indice SSN et Aa avec asymptote en 3π 2 . . . . . . 11 Indice Aa lissé en fonction du SSN lissé. . . . . . . . . . . . . . . . . . . . . . 12 Coefficient d'optimisation a et b entre 1868 et 2010 et coefficient de corrélation 13 Graphique du modèle optimisé tous les 12 mois . . . . . . . . . . . . . . . . . 14 Résultats de l'optimisation de la tendance des coefficients a et b pour le Modèle Aa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Modèle avec les valeurs de coefficients estimées entre 1868 et 2010 et extrapolées pour la période entre 2010 et 2019. . . . . . . . . . . . . . . . . . . . . . 16 Résidus du modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Ajout des dernières données de Aa mensuel (janvier 2010 -juin 2011). . . . . . 18 Nombre d'occurrences annuelles et mensuelles des orages magnétiques. . . . . 19 Distribution de l'indice Dst (log-transformé). . . . . . . . . . . . . . . . . . . 20 Représentation du F10.7 et Dst réunis . . . . . . . . . . . . . . . . . . . . . . 21 Graphique de la densité de probabilité d(Dst ≤ -100nT ) en fonction du F10.7 moyen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Liste des tableaux 1 Bilan des différentes sources d'erreurs GPS . . . . . . . . . . . . . . . . . . . Ordre de grandeur de chaque indice . . . . . . . . . . . . . . . . . . . . . . . Bilan des coefficients de corrélation entre les différents indices . . . . . . . . .

2. 2

 2 Etat de l'art du GNSS (N. Suard -Cours Supelec) Pour comprendre comme il se doit l'objet de l'étude décrit ici, il est essentiel de préciser les différents principes de navigation et propagation dans le domaine Spatial. Le GPS (Global Positioning System) a été expérimenté dès 1978 mais l'ensemble des satellites mis en orbite n'ont été déclarés en service qu'à partir de 1995. D'autres systèmes accompagnent le GPS tel que GLONASS (Système de positionnement par satellite Russe), GALILEO (Européen) le plus récent.Il existe plusieurs familles d'utilisateurs de ces systèmes de positionnement, les services GPS répandus tel que délivrés au travers de récepteurs applicatifs type Tomtom ou Garmin sont pour le grand public (civil) ; ce service est le moins complexe à mettre en oeuvre et demande le moins de ressources du fait que la grande majorité des utilisateurs n'évolue que sur une composante pseudo-horizontale. Ensuite nous avons les applications militaires de ces systèmes, aujourd'hui il existe une forte dépendance aux USA des différents pays développés qui ne possèdent pas encore de système de navigation par satellite totalement indépendant. Chacun cherche à devenir indépendant des autres dans l'éventualité d'un conflit majeur. Le monde scientifique utilise beaucoup ces services pour toutes les études spatiales en cours et futures. D'autres utilisations de ces systèmes sont déjà des niches économiques, nous avons l'exemple des agriculteurs qui utilisent de plus en plus ces services afin d'obtenir les meilleurs rendements possibles.

  et il serait couplé à un cycle de 22 ans. Selon les études précédentes(Rangarajian, 1997 & Clúa de Gonzalez, 2000) l'analyse des indices par intervalles permet l'étude d'occurrence des extrémums. Echer, 2011 a réalisé une étude statistique sur les occurrences des orages magnétiques sévères à l'aide d'une distribution exponentielle de l'indice Dst.[START_REF] Srivastava | A logistic regression model for predicting the occurrence of intense geomagnetic storms[END_REF][START_REF] Srivastava | The Challenge of Predicting the Occurrence of Intense Storms[END_REF] a mis en place un modèle de régression logistique qui peut par apprentissage, prédire l'occurrence des orages sévères et intenses à l'aide de nombreuses variables. Le modèle de Cox a été expérimenté sur l'occurrence des orages de fortes intensité, c'est un modèle basé sur les hasards proportionnels[START_REF] Azaïs | Etudes des orages magnétiques, Méthode d'évènements rares pour l'intégrité et la continuité EGNOS -GALILEO[END_REF], ils en arrivent à la conclusion que l'intensité d'apparition d'orages magnétiques dépend essentiellement de la position dans le cycle, l'intensité atteint son maximum en phase descendante de cycle indépendamment du niveau d'orage magnétique.2.5 Objectifs du stageAu vu des informations précédentes, il m'a été demandé à partir des données décrites ci après de valider et quantifier autant que peut se faire les différentes relations bibliographiques exposées et éventuellement d'en mettre en évidence de nouvelles, ceci à l'aide d'outils et de méthodes à développer/proposer par moi-même.Les données de mesures de l'activité solaire sont recensées régulièrement depuis les années 1750, et c'est seulement à partir des années 1840 qu'un amateur d'astronomie Heinrich Schwabe, qui cherchait à repérer une éventuelle planète dans l'orbite de Mercure, commence à compter de manière systématique les taches du soleil. Il espérait apercevoir cette planète en la différenciant des autres taches. C'est en réalisant ces comptages qu'il évoqua la forte probabilité que l'occurrence de ces taches présente un cycle moyen de 11ans. Pendant la même période Rudolf Wolf s'intéressa également aux travaux de H. Schwabe et améliora la qualité de comptage par estimation du nombre de taches[START_REF] Hathaway | The shape of the sunspot cycle[END_REF]. Dans les années 1900, Edward Maunder s'intéressa à la période particulièrement froide entre 1650 et 1720 et mit en évidence une très faible activité solaire durant cette période[START_REF] Mendoza | Geomagnetic activity and wind velocity during the Maunder minimum[END_REF]. 1958, Gleissberg met en évidence un cycle de 80ans dans les taches solaires. Les taches solaires ont fait l'objet de nombreuses observations durant les siècles derniers. C'est pour cela que nous avons une très grande période de relevés de cet indice solaire. Depuis que les données régulières de l'activité solaire sont accessibles nous en sommes au 24ième cycle.Nous avons également un autre indice traduisant l'activité solaire, le F10.7 : il mesure le flux radio à 10.7 cm de longueur d'onde, ce flux provient de la région chromosphérique du soleil et dans la basse couronne suite à l'occurrence de taches solaires. Les niveaux d'intensité radio sont mesurés en trois points de l'héliosphère :-Sur une région non-perturbée de la surface, -Sur une région active en développement, -Sur une région intense et plus forte que le niveau d'activité journalier.

  FIGURE 1 -Ensemble des indices solaires et géomagnétiques en fonction du temps.

  Densité Spectrale des Puissances (DSP) La DSP est couplée avec la transformée de Fourier sur les signaux, elle présente le spectre du signal avec en abscisse sa fréquence et en ordonnée les puissances de chaque fréquence. La transformation est réalisée à l'aide d'une fonction d'autocorrélation (équation 3 pour tout t) et d'une transformée de Fourier sur cette même fonction (équation 4). Ce spectre nous permet de déceler une plage de fréquence locale qui présente une forte puissance. Il est alors intéressant de sélectionner la plage qui nous intéresse, c'est-à-dire mettre à zéro les fréquences qui ne nous intéressent pas. Après cette étape, nous pouvons recomposer le signal et voir quel phénomène est mis en évidence. La DSP est un outil indispensable à l'analyse des spectres des différentes variables.

  par ondelettes La méthode par ondelettes est utilisée depuis les années 1930. Les ondelettes font partie de ces méthodes spécifiques aux signaux non-stationnaires qui forment une grande majorité comparés aux signaux stationnaires. Ces méthodes sont divisibles en deux ensembles, les méthodes temps-fréquence : avec une prise en compte d'une évolution temporelle du contenu fréquentiel d'un signal ; et les méthodes temps-échelle (ondelettes) : qui décrivent l'évolution temporelle d'un signal à une échelle d'observation. (P. Flandrin et Y. Meyer, Le courrier du CNRS n˚77). Les ondelettes sont utilisées pour leur facilité d'utilisation et d'interprétation. Ses résultats dépasseraient l'efficacité de la transformée de Fourier. Cependant il ne faut pas voir la transformée de Fourier et Ondelette comme deux méthodes distinctes mais bien complémentaires. Par exemple, la transformée de Fourier apparaît comme un prétraitement pour isoler un son, l'ondelette intervient en suivant pour affiner l'isolement.La transformée par ondelettes bien que particulière présente la même base Hilbertienne que la transformée de Fourier. L'ondelette « mère » ou de base que l'on appliquera sur plusieurs échelles selon les coefficients a et b se présente sous la forme de l'équation 5. La transformée est présentée en équation 6 et l'équation 7 est la représentation de la décomposition des signaux.

  Transformée en ondelette continue (CWT) et Image processing La transformée en ondelette continue nous permet de travailler sur n'importe quelle échelle, de cette manière nous pouvons produire un gradient de courbe selon l'échelle utilisée. Tous les vecteurs de coefficients créés grâce à cette transformée pourront être réutilisés par un outil de Scilab très utile d'un point de vue descriptif, c'est une boîte à outils du logiciel : Scilab Image Processing (SIP Toolbox). Elle permet de réutiliser l'ensemble des courbes dont nous parlions précédemment pour les mettre sous formes d'images avec un gradient de couleurs.Nous pouvons plus aisément comparer les coefficients de transformation entre les indices et la manière dont réagit l'ondelette mère aux signaux. Par cette méthode nous pouvons observer facilement si deux signaux suivent un processus commun à différentes échelles.Transformée en ondelette discrète (DWT)C'est par abus de langage qu'elle est appelée de cette manière, c'est en fait la forme discrétisée de la transformée en ondelette continue. Nous l'utilisons car nous aurons besoin de recomposer les signaux après traitement par les ondelettes. Le choix de la DWT s'est fait sur elle car elle effectue un traitement équivalent à la CWT et l'inverse de la DWT (IDWT) est disponible sur Scilab. L'inverse de la CWT est impossible due à la redondance des coefficients dans la décomposition sur toutes les échelles. Cependant cette DWT présente deux contraintes, elle ne peut traiter les données que sur des intervalles d'échelles dyadique (par exemple : 264-1024), de plus elle ne peut pas prendre en charge les ondelettes continues telle que l'ondelette de Morlet.

L

  FIGURE 2 -(a) Traitement de l'indice Aa par la transformée de Fourier(gauche) et (b) DSP du signal avec la fenêtre sélectionnée (droite).

4. 3 FIGURE 5 -

 35 FIGURE 5 -Graphiques représentant les indices SSN et Aa mensuels en fonction du temps.

FIGURE 6 -FIGURE 7 -

 67 FIGURE 6 -Décomposition multi-échelle par l'ondelette « db5 » sous forme de graphique (à gauche) et sous forme d'image (à droite), 1024 points en 256 échelles.

FIGURE 8 -

 8 FIGURE 8 -Graphiques de corrélation entre SSN et Aa sur l'ensemble des échelles après transformation par ondelette de Morlet à gauche (a), et Daubechies 5 à droite (b).

4. 4 FIGURE 9 -

 49 FIGURE 9 -Graphique des signaux recomposés après ondelette de Daubechies 5 entre l'échelle 64 et 1024, R=0,72 pour la courbe rouge.

FIGURE 16 -

 16 FIGURE 16 -Analyse des résidus du modèle Aa par rapport aux observations.

FIGURE 17 -

 17 FIGURE 17 -Ajout des dernières données de Aa mensuel (janvier 2010 -juin 2011).

FIGURE 20 -

 20 FIGURE 20 -Représentation du F10.7 et Dst réunis, les barres verticales représentent les sections que nous avons choisis.

FIGURE 21 -

 21 FIGURE 21 -Graphique de la densité de probabilité d(Dst ≤ -100nT ) en fonction du F10.7 moyen centré par la moyenne générale (127.5 SFU).

  Cette relation a été mise en évidence grâce à plusieurs techniques de traitement de signal, la transformée de Fourier nous a permis de connaître les corrélations existantes entre les différents indices, les transformations multi-échelles par les ondelettes ont mis en évidence les ressemblances des signaux solaires et géomagnétiques. Nous avons défini ensuite une équation pouvant relier les deux signaux lissés (SSN et Aa) entre eux, c'est une équation linéaire dont les deux paramètres varient au cours du temps selon des sinusoïdes connues. Différents travaux ont démontré une relation également linéaire entre activité solaire et géomagnétique (Ortiz de

  présentent un décalage de 1an sur les signaux pour réaliser l'analyse. De ce fait nous avons un moyen de prédiction à moyen-terme de l'activité géomagnétique sur l'indice Aa. L'analyse des résidus montre leur normalité et leur homoscédasticité ce qui valide le modèle. Grâce à ce modèle et le décalage existant entre les deux indices nous avons la possibilité de prévoir l'activité géomagnétique pour les 16 mois à venir et même pour la période du cycle solaire 25 (2008-2020) grâce aux données du MSFC. • Les occurrences des forts évènements Beaucoup de travaux ont été effectués sur l'occurrence des forts orages magnétiques afin de les prévoir, cet objectif est très important par rapport à la finalité de l'étude qui est la sauvegarde des données GPS. Azaïs et al., 2011 ont établi l'intensité d'apparition de plusieurs niveaux d'orages magnétiques, elle dépend principalement de la position dans le cycle solaire. Comme dans leurs résultats, nous voyons plus d'occurrences sur la seconde partie des cycles. Nous avons pu définir la probabilité d'occurrence à posteriori des orages dit « sévères » dont le Dst était inférieur à -100nT, elle est d'environ 2,24% sur les 50 dernières années. La distribution de l'indice Dst est de type log-normale ce qui permet d'affirmer que nous aurons toujours ce risque dans le futur fut-il aussi infime. Echer et al., 2011 ont employé la même méthode dans le but de déterminer cette probabilité d'occurrence, en utilisant une distribution exponentielle. Nous avons mené la même analyse sur chaque cycle de ≈ 1965 à 2011, et nous révélons une relation linéaire très significative entre la probabilité d'occurrence d'orages magnétique (Dst ≤ -100nT ) et l'activité moyenne des cycles solaires. Nous avons observé les apparitions d'orages par rapport aux cycles solaires et ceux-ci apparaissent préférentiellement dans la seconde partie des cycles. La prédiction de tels évènements est une tâche plus longue et demande un travail très approfondi sur le type de technique que nous devrions utiliser. L'utilisation des méthodes d'apprentissage tel que le Random Forest ou les Réseaux de Neurones peuvent être très efficaces, mais les variables explicatives de l'activité solaire ne permettent pas aux modèles d'apprendre assez pour qu'ils puissent ensuite prédire l'occurrence des orages. Dans notre cas ces méthodes peuvent être très puissantes, mais nous sommes limités par la qualité des variables d'entrées qui ne permettent pas d'effectuer une prédiction acceptable. Srivastava, 2005 et 2006 ont défini un modèle de prédiction par une simple régression logistique, les variables d'entrées sont nombreuses, elles sont séparées en deux catégories : les variables solaires avec plusieurs composantes des CMEs et les variables interplanétaires telles que le champ magnétique interplanétaire et les vents solaires (densité et vitesse). La plupart des travaux réalisés sur ce sujet, attribuaient à chaque orage un type d'évènement solaire, il peut y avoir une dizaine d'évènements différents qui sont la cause des orages magnétiques.

  angulation (position x, y, z sur un repère centré Terre) et un dernier pour corriger l'erreur de rattachement au temps système. C'est grâce au satellite SBAS (Satellite Based Augmentation System) que l'intégrité (une des quatre performances GNSS que nous développerons plus tard) des données est assurée. Or il arrive que les informations GPS soient dégradées régulièrement et donc l'intégrité de la position GPS est défaillante. Il est fort probable que ces problèmes soient dus aux évènements géomagnétiques. Ces évènements géomagnétiques sont très souvent les conséquences d'une forte activité solaire, mais cette affirmation n'est pas prouvée statistiquement.

	2 Contexte, problématique
	2.1 Présentation de l'organisme d'accueil
	Mon stage s'est déroulé au Centre National des Etudes Spatiales (CNES Toulouse) au sein
	du service Localisation-Navigation (LN) et Instrumentation, Télémesure/Télécommande, Pro-
	pagation (ITP), ces services font eux même partie de la sous-direction Radio-Fréquence (RF).
	J'ai pu travailler en collaboration avec Norbert Suard (Ingénieur, LN) et Françoise Carvalho
	(Ingénieur, ITP).
	Le CNES est un établissement public à caractère industriel et commercial (EPIC) placé sous
	la double tutelle du ministère de l'Enseignement et de la Recherche et du ministère de la Dé-

L'ionosphère est une partie de l'atmosphère terrestre située entre 60 et 1000 à 1200 km d'altitudes. C'est le lieu où de nombreuses réactions photochimiques dues aux radiations solaires se produisent. Celles-ci agissent sur les gaz présents ce qui produit le plasma ionisé. Ces ions sont alors piégés dans le champ magnétique de la terre. Plusieurs évènements magnétiques peuvent se produire suite à ces réactions dont les plus connus sont les Aurores boréales au niveau du pôle Nord. Cela peut aussi induire des orages magnétiques plus ou moins sévères, les orages de forte puissance peuvent entraîner un effet de Plume qui est la conséquence de particules chargées créant une variation dans le champs magnétique assez forte pour se propager suivant les lignes magnétiques, cela apparaît sur les cartes sous forme de plume magnétique d'où son nom (Eurocontrol's funded study -2010). La relation entre les indices magnétiques et les évènements ionosphériques semble être évidente (Kp et occurrence d'aurores boréales corrélés), or dans cette étude nous voulons mettre en évidence les relations qui peuvent exister entre l'activité solaire et l'activité géomagnétique. Les évènements que nous avons énumérés précédemment semblent prouver cette relation mais nous ne pouvons pas l'affirmer statistiquement et de manière significative pour le moment. Si tel était le cas et en fonction des prévisions que l'on pourrait faire sur l'activité solaire, nous pourrions en déduire des prévisions sur les variations des indices magnétiques. Dans un premier temps nous ferons le descriptif de l'organisme d'accueil, de la problématique générale ainsi que de l'ensemble des données dont nous disposons. Dans une seconde partie nous ferons la liste exhaustive des méthodes utilisées pour atteindre les objectifs de notre étude. Dans la partie suivante nous détaillerons les résultats que nous avons obtenus et nous terminerons par une discussion des résultats et une conclusion. fense. Fondé en 1961 son siège social est situé à Paris ; il existe trois centres spatiaux : le centre spatial d'Evry appelé la direction des lanceurs, celui de Toulouse spécialisé dans les véhicules et systèmes orbitaux, et le centre en Guyane française, dit le port spatial de l'Europe. A l'origine, le CNES avait pour mission de proposer une stratégie spatiale et la mettre en oeuvre. Il était capital à cette époque d'après guerre de se mettre au même niveau que les Etats-Unis et l'URSS pour ce qui est des avancées spatiales. En association avec l'ESA (Agence Spatiale Européenne), le CNES a permis de mettre en place des projets de grandes envergures, le premier élément majeur fut le développement du lanceur Européen : Diamant A dont le premier lancement a été réalisé le 26 Novembre 1965 pour mettre en orbite le premier satellite Européen Astérix. Cet événement mit la France à la troisième position des puissances mondiales dans le domaine spatial. Dès lors, le CNES n'a pas cessé de développer les instruments d'exploration spatiale et de navigation.

  'activité solaire afin d'avoir une meilleure précision.Ahluwalia, 2000 a prédit un maximum de SSN à 131,5 pour le cycle 23, or ce cycle a présenté un SSN maximum moyenne de 120. C'est ici que nous pouvons remarquer la difficulté d'un tel objectif. Le cycle solaire ne semble pas être régulier malgré sa forme périodique, nous pouvons nous baser sur certaines observations fondées pour prédire une activité plus faible pour le début de notre siècle, en effet Silverman en 1992 affirme qu'à chaque début de siècle nous pouvons observer une activité solaire minimale et prolongée. Un modèle qui respecte bien la forme des cycles précédents a été déterminé[START_REF] Hathaway | The shape of the sunspot cycle[END_REF] en utilisant seulement 4 paramètres, l'amplitude du cycle suivant est estimée grâce à la longueur du cycle précédent, une des difficultés de ce modèle est la détermination du t0 de chaque cycle. Quelques études sur la prédiction du prochain cycle font l'état des différents travaux effectués et donnent le résultat sur

	E. Rifa	2 CONTEXTE, PROBLÉMATIQUE
	2.4 Recherche bibliographique besoin de mieux quantifier l
	Avec une base de documents fournis au démarrage du stage, j'ai pu prendre connaissance
	du sujet en profondeur, des notions fondamentales de la navigation jusqu'à la prédiction des
	prochains cycles solaires. Les recherches principales que j'ai menées ont été fructueuses et réa-
	lisées sur les bases de données telles que ScienceDirect ou même Pascal qui fournissent des
	publications scientifiques dans tous les domaines dont celui du spatial.
	2.4.1 Relations
	Plusieurs études ont été réalisées afin d'obtenir une relation entre l'activité solaire et un ou
	plusieurs indices terrestres. Il est établi qu'il y a une relation linéaire entre les indices ionosphé-
	riques et l'activité solaire (Ortiz de Adler, 1997) et il existe une relation entre la concentration
	en électron de la région D de l'ionosphère et l'activité solaire (Danilov, 1998), ces études per-
	mettaient de donner un premier lien entre l'activité solaire et les différentes activités ionosphé-
	riques. Une étude publiée en 2000 par H.S. Ahluwalia décrit une relation entre le nombre de
	taches solaires et un indice géomagnétique (Ap), ils ont pu définir une relation entre ces deux
	indices avec une corrélation proche de 0,82 avec les données annuelles. Très récemment une
	relation entre l'activité solaire et un paramètre ionosphérique : le Fof2 a été mise en évidence et
	celle-ci semble être changeante, entre 1975 et 2000 elle était qualifiée de normale et à partir de
	2001-2002 il y eut des changements dans les coefficients (Lukianova, 2011).
	Une étude très proche de notre problématique a été réalisée par G.K. Rangarajian en 1997
	Comme nous l'avons expliqué dans le paragraphe précédent, ces systèmes de correction et où il trouve une corrélation très significative entre les Kp=4-5 et l'activité solaire. L'analyse a
	le GPS seuls présentent une grande sensibilité aux évènements géomagnétiques. C'est dans ce été faite sur l'occurrence des différents intervalles de Kp. Dans le même objectif une relation
	cadre que s'insère l'étude des relations entre les indices solaires et géomagnétiques (sujet du entre l'indice Aa et les taches solaires a été analysée (Clùa de Gonzalez et al., 2000). De la
	stage) puisqu'il y a un très grand intérêt à connaître cette relation. Le cas idéal serait de pouvoir même manière l'indice Aa a été étudié par intervalles.
	estimer l'activité géomagnétique, qui est traduite à travers la mesure des différents indices,
	seulement à partir de la mesure de l'activité solaire. 2.4.2 Prédictions
	Une fois cette relation établie nous chercherons à définir l'activité géomagnétique future et Les scientifiques ont publié de nombreux articles sur le sujet et la plupart cherchent depuis
	comme autre objectif : trouver la probabilité d'occurrence d'orages géomagnétiques sévères. une vingtaine d'années à définir le modèle qui permettrait de prédire la forme des cycles so-
	laires. Nous savons d'ors et déjà que le cycle solaire n'a pas une forme symétrique, les phases
	ascendantes et descendantes diffèrent par leur pente et il semblerait qu'à chaque fin de cycle,
	un processus de randomisation se mette en place (Kim, 2006) rendant difficile la prédiction du
	cycle solaire suivant à partir du précédent. Les publications du début des années 2000 parlent
	d'un cycle 24 (début en 2008) beaucoup plus faible que le cycle 23 (Schatten, 2003) mais il y a

le maximum de taches solaires à venir

[START_REF] Echer | Prediction of solar activity on the basis of spectral characteristics of sunspot number[END_REF][START_REF] Gholipour | Solar activity forecast : Spectral analysis and neurofuzzy prediction[END_REF]

; Rigozo et al., 2010), l'ensemble des études dont Echer et al., 2004 tendaient à sur-estimer le maximum du cycle 23.

TABLE 2 -

 2 Pour définir la relation indice solaire-géomagnétique, nous avons fait la sélection de plus de 19000 données journalières par indice ce qui nous permet de faire l'analyse entre 1957 et 2010 sur tous les indices disponibles. En effet les données de Dst ne sont disponibles que Ordre de grandeur de chaque indice Nous avons un jeu de données discontinu qui ne nous permet pas d'utiliser les méthodes de traitement que nous voulons mettre en oeuvre. Celles-ci demandent un signal continu sans données manquantes. Afin de pallier ce problème nous avons réalisé une interpolation par moyenne mobile, qui aura moins d'impact sur les données qu'une interpolation linéaire. La moyenne mobile (équation 1) tiendra compte de la tendance locale alors qu'une simple interpolation linéaire relie deux points éloignés l'un de l'autre par une droite.Sur les données brutes nous pouvons deviner certaines tendances qui peuvent être communes entre indices solaires et géomagnétiques, il est alors important de définir une stratégie afin d'observer d'éventuelles liaisons temporelles entre ces variables. L'idée serait d'analyser les données dans un sous-espace afin de déterminer les échelles des signaux bruts qui présentent des processus commun entre les différentes variables.La transformée de Fourier permet facilement de trouver les variations régulières de signaux stationnaires. Malgré la non-stationnarité de nos signaux, notamment celui des taches solaires, nous pouvons récupérer des signaux de hautes puissances qui permettraient de mettre en évidence certains cycles.

	depuis janvier 1957. Nous avons la possibilité d'effectuer l'analyse sur les données journalières,
	mensuelles et annuelles.
	Nous avons choisis d'analyser les signaux SSN et F10.7 sous leur forme brute, tels qu'ils
	sont mesurées, de cette manière nous ne nous éloignerons pas des valeurs réelles ; les in-
	dices géomagnétiques sont moyennés journalièrement pour l'analyse. Nous avons réalisé une
	moyenne mobile (équation 1) sur tous les signaux pour une première approche par analyse des
	Indices Minimum Maximum
	SSN	8	355
	F10.7	64,5	461,8
	Kp	0	83,4
	Ap	0	280
	Aa	0	352,1
	Cp	0	9
	C9	0	3,3
	Dst	-35	240
	Le nombre de taches solaires correspond au nombre observé par jour, or les indices tels que
	le Kp ou le Ap sont mesurés toutes les 3 heures. De cette manière une étude peut également se
	faire sur les Kp moyen, maximum et peut être minimum par jour. Une étude précédente met en

corrélations entre les indices. La moyenne mobile a été appliquée notamment pour l'indice SSN qui présente de nombreuses plages de données manquantes, l'indice F10.7 et Aa présentaient eux aussi quelques valeurs manquantes. Le Dst étant majoritairement négatif nous l'avons rendu positif, pour pouvoir le comparer. Le tableau 2 résume les ordres de grandeur des différents indices.

avant la relation entre les Ap et les taches solaires avec un décalage de un an

[START_REF] Ahluwalia | Solar cycle 23 prediction update[END_REF]

. Pour notre part, nous avons également observé un tel décalage entre les indices géomagnétiques et l'indice solaire mais de plus d'un an (480 jours en moyenne, 16 mois). Ce décalage moyen a été estimé par l'observation des pics d'indices moyens et par une analyse de phase. [section 4.4] avec i : le ième élément, k : longueur de la demi-fenêtre.

3.3 Généralités sur le traitement du signal L'observation des phénomènes naturels nécessite des méthodes efficaces qui permettent de minimiser les erreurs de mesures ainsi que de préserver l'intégrité des informations. Mais l'interprétation des observations ne se fait pas sans un traitement au préalable, car même avec la précision la plus parfaite, il existe toujours une erreur standard, un bruit. Il ne faut pas oublier que nous travaillons sur des données qui mesurent des phénomènes naturels et par définition la nature représente le chaos. Cependant, l'environnement est composé de différents compartiments qui peuvent interagir entre eux et des relations peuvent être mises en évidence. Les principales méthodes que nous utilisons pour cette étude proviennent d'une matière scientifique communément appelée : le traitement du signal. « C'est l'ensemble des méthodes qui permet de transformer les données dans un but final d'exploitation ; c'est l'extraction du maximum d'informations utile d'un signal bruité. » (Courrier du CNRS, n˚77). Ce domaine d'étude est très riche en ressources ainsi qu'en domaines d'applications. Dans chaque relevé de données, nous avons une abondance d'information qui ne nous permet pas de tirer des conclu-Dans un premier temps nous utilisons la Fast Fourrier Transform (FFT) ou la transformée de Fourier rapide (équation 2). Elle permet de convertir les signaux dans le domaine des exponentielles complexes sur une base Hilbertienne ou base orthonormale (BON).

  La spline est une méthode de lissage de signal plus directe que celle par l'utilisation du traitement du signal. Avec cette technique nous pouvons obtenir n'importe quel type de lissage des signaux, il suffit de choisir le degré de liberté qu'il faut. Mais contrairement à la DWT nous ne sommes pas capables de savoir à quelle échelle nous nous situons sur le signal lissé. Si nous comparons les signaux reconstruits après ondelette avec les résultats de spline et que nous 'avantage d'une telle méthode est que nous pouvons introduire autant de paramètres que l'on veut pour définir le meilleur modèle. Elle estime les paramètres afin d'obtenir le meilleur ajustement, celui-ci se réalise par la méthode des moindres carrés. Nous définissons alors une fonction f à n paramètres résumés en un vecteur ω , nous cherchons une relations entre x, y qui ont chacun une longueur k : y = f (x, ω) , et nous ajustons les valeurs du vecteur ω tel que :

		k
	SCE =	i=1 ∑ (y i -f (x i , ω)) 2 →≈ 0	(8)
	Nous intégrons des coefficients à ajuster proche des vrais coefficients afin que la fonction
	d'optimisation converge vers le bon minimum local. C'est l'inconvénient de cette technique,
	obtenons un bon coefficient de corrélation, nous pourrons estimer avoir le bon signal. Nous
	choisissons d'utiliser cette technique car nous n'avons pas les inconvénients des effets de fin de
	signal avec la spline. En effet la DWT a tendance à fausser la fin du signal alors que la spline

suit les données jusqu'à la dernière. La méthode des splines est disponible sous le logiciel R avec la librairie « splines ». 3.6 Optimisation, régression non-linéaire Lselon la fonction et les coefficients que nous intégrons, il est possible que le gradient d'optimisation présente plusieurs minimums d'écarts au carré.

TABLE 3 -

 3 Bilan des coefficients de corrélation entre les différents indices

	Fr 0,0001-0,006 SSN F10.7 Kp	Ap	Aa	Cp	C9	Dst
	SSN	1	0,98 0,28 0,31 0,25 0,27 0,26 0,34
	F10.7	0,98	1	0,28 0,32 0,27 0,27 0,26 0,36
	Kp	0,28 0,28	1	0,94 0,94	1	1	0,72
	Ap	0,31 0,32 0,94	1	0,95 0,94 0,93 0,78
	Aa	0,25 0,27 0,94 0,95	1	0,94 0,94 0,78
	Cp	0,27 0,27	1	0,94 0,94	1	1	0,72
	C9	0,26 0,26	1	0,93 0,94	1	1	0,72
	Dst	0,34 0,36 0,72 0,78 0,78 0,72 0,71	1
	nous l'attendions. Ensuite les indices géomagnétiques sont en grande majorité bien corrélés à
	l'exception du Dst qui présente moins de 80% de corrélation avec les autres indices géoma-
	gnétiques. Une étude a été réalisée sur les différences entre les indices Ap et Aa (Coffey et al.,
	2001), la plupart des décorrélations ont été décelées par l'analyse des indices par intervalles et
	par horaires, mais de manière générale nous avons une corrélation approchante à celle que nous
	trouvons ici.

Afin de réaliser l'étude la plus efficace et complète nous choisissons de n'utiliser que l'indice Aa des indices géomagnétiques. Cet indice présente le meilleur compromis et la plus large période de recensement (entre 1868 et 2011) contrairement aux autres indices qui ne sont disponibles que depuis 1932 ou 1957 pour le Dst.

  La figure 4 met en évidence les différents pics de puissances sur la densité spectrale des puissances. Le pic le plus puissant que nous observons est celui du cycle de 11 ans sur les deux signaux. Celui-ci est simplement dû à l'activité solaire, nous le retrouvons bien dans le géomagnétisme. Ensuite nous observons un cycle de 5,5 ans, celui-ci serait une harmonique du cycle de 11 ans(Prestes et al., 2005), cependant nous pouvons émettre l'hypothèse que son origine viendrait du changement dans la dynamo solaire. Nous observons aussi un pic à 2,7 ans, remarqué aussi parPrestes et al., 2005 ; ce processus est dû à une structure à double ou triple pics pendant la phase descendante du cycle solaire. Les deux pics suivant à environ 22 et 17 mois ressortent du géomagnétisme seulement, seul le pic de 17 mois est décelé aussi dans l'étude dePrestes et al., 2005 mais ne donne pas de précisions. Ensuite un pic à 6 mois est mis en évidence sur la figure4qui n'apparaît encore que sur le géomagnétisme, ce cycle est simplement dû au rythme saisonnier de la Terre. Nous voyons maintenant un pic à 27 jours sur les signaux solaires et géomagnétique, ce cycle a été mis en évidence par[START_REF] Bartells | Twenty-seven day recurrences in terrestrial-magnetic and solar activity, 1923-1933[END_REF] ; le temps de ré-

). Ces sauts de puissances signifient qu'un processus régulier est présent tout le long du signal. FIGURE 4 -Spectre des Puissances des signaux SSN (haut) et Aa (bas).

volution du soleil de 27 jours environ est à l'origine de ce cycle. Enfin nous pouvons voir une forte puissance du pic de 13 jours environ, ce pic ne ressort pas autant sur les signaux solaires, nous n'avons pas d'information sur l'origine de ce pic de puissance.

Le jeu de données actuel ne nous permet pas d'observer l'ensemble des pics avec la méthode que nous utilisons. Grâce à la méthode du Multitaper

(Prestes et al., 2005)

, ils ont pu déterminer certains cycles que nous n'observons pas tel que le cycle de 80-90 ans sur l'activité solaire, et une période d'environ 35 ans sur l'activité géomagnétique qui n'apparaît pas dans l'activité solaire.

Le cycle de Gleissberg (80-90 ans) n'apparaît pas sur notre spectre de puissance car les données dont nous disposons ne permettent pas de mettre en évidence celui-ci. Střeštík, 2009 démontre à l'aide d'une analyse spectrale la présence d'une périodicité de 1,3 ans visible sur l'indice Ap et non sur le SSN entre 1932 et 2007. Cette période n'apparaît comme un processus visible que sur un laps de temps de l'activité solaire.

Remerciements


Adrien Le Coënt
email: adrien.le-coent@ens-cachan.frdevuyst@cmla.ens-cachan.fr

Laurent Fribourg
email: fribourg@lsv.ens-cachan.fr

Nicolas Markey
email: markey@lsv.ens-cachan.fr

Florian De Vuyst

Ludovic Chamoin
email: chamoin@lmt.ens-cachan.fr

Distributed Synthesis of State-Dependent Switching Control

We present a correct-by-design method of state-dependent control synthesis for linear discrete-time switching systems. Given an objective region R of the state space, the method builds a capture set S and a control which steers any element of S into R. The method works by iterated backward reachability from R. More precisely, S is given as a parametric extension of R, and the maximum value of the parameter is solved by linear programming. The method can also be used to synthesize a stability control which maintains indefinitely within R all the states starting at R. We explain how the synthesis method can be performed in a distributed manner. The method has been implemented and successfully applied to the synthesis of a distributed control of a concrete floor heating system with 11 rooms and 2 11 = 2048 switching modes.

Introduction

The importance of switched systems has grown up considerably these last years in reason of their ease of implementation for controlling cyber-physical systems. A switched system is a family of sub-systems, each with its own dynamics characterized by a parameter mode u whose values are in a finite set U (see [START_REF] Liberzon | Switching in systems and control[END_REF]). However, due to the composition of many switched systems together, the global switched system has a number of modes and dynamics which increases exponentially. Take for example a heating system for a building of 11 rooms (see [START_REF] Kim G Larsen | Online and Compositional Learning of Controllers with Application to Floor Heating[END_REF]): each room i has a heater with 2 mode values {off,on}. This makes a combination of 2 11 = 2048 mode values. If we want to analyze the evolution of a trajectory on a horizon of K units of discrete time, we have to consider the dynamics corresponding to 2 11K possible sequences of modes, which is intractable even for small values of K. It is therefore essential to design compositional methods in order to obtain control methods of switched systems that give formal guarantees on the correct behavior of the cyber physical systems.

In this paper, we give a symbolic compositional method which allows to synthesize a control of linear discrete-time switched systems that is guaranteed to satisfy attainability and stability properties.

The method starts from an objective region R of the state space, which is rectangular (i.e., is a product of closed intervals of reals). It then generates in a backward manner, using linear programming techniques, an increasing sequence of nested rectangles {R (i) } i≥0 such that every trajectory issued from R (i) is guaranteed to reach R (i-1) in a bounded number of time units. Once R (0) = R is reached, the trajectory is also guaranteed to stay in R indefinitely (stability). The method relies on a simple operation of tiling of the rectangles R (i) in a finite number of sub-rectangles (tiles), using a standard operation of bisection. Although the method works in a backward fashion, it does not require to inverse the linear dynamics of the system (via matrix inversion), and does not compute predecessors of symbolic states (tiles), but only successors using the forward dynamics. This is useful in order to avoid numerical imprecisions, especially when the dynamics are contractive, which happens often in practical systems (see [10]).

Another contribution of this paper is a technique of state over-approximation which allows a distributed control synthesis: this over-approximation allows subsystem 1 to infer a correct value for its next local mode u 1 without knowing the exact value of the state of sub-system 2. This distributed synthesis method is computationally efficient, and works in presence of partial observability. This is at the cost of the performance of the control which usually makes the trajectories attain the objective area in more steps than with a centralized approach.

Related Work

In symbolic analysis and control synthesis methods for hybrid systems, the method of backward reachability and the use of polyhedral symbolic states, as used here, is classical (see, e.g., [START_REF] Asarin | Effective synthesis of switching controllers for linear systems[END_REF][START_REF] Jeremy H Gillula | Applications of hybrid reachability analysis to robotic aerial vehicles[END_REF]). The use of tiling or partitioning the state-space using bisection is also classical (see, e.g., [START_REF] Jaulin | Applied Interval Analysis with Examples in Parameter and State Estimation[END_REF][START_REF] Girard | Low-complexity switching controllers for safety using symbolic models[END_REF]). The main original contribution of this paper is to give a simple technique of over-approximation, which allows one component to estimate the symbolic state of the other component, in presence of partial information. This is similar in spirit to an assumeguarantee reasoning where the controller synthesis for each sub-systems assumes that some safety properties are are satisfied by the others [START_REF] Alur | Reactive modules[END_REF]9]. In contrast to [START_REF] Fribourg | Game-based Synthesis of Distributed Controllers for Sampled Switched Systems[END_REF], we do not need, for the mode selection of a sub-system, to explore blindly all the possible mode choices made by the other sub-system. This yields a drastic reduction of the complexity 4 . This approach allows us to treat a real case study which is intractable with a centralized approach. This case study comes from [START_REF] Kim G Larsen | Online and Compositional Learning of Controllers with Application to Floor Heating[END_REF], and we use the same decomposition of the system in two parts (rooms 1-5 and rooms 6-11). In contrast to the work of [START_REF] Kim G Larsen | Online and Compositional Learning of Controllers with Application to Floor Heating[END_REF] which uses an on-line and heuris-tic approach with no formal guarantees, we use here an off-line formal method which guarantees attainability and stability properties.

Implementation

The methods of control synthesis both in the centralized context and in the distributed context have been implemented in a prototype written in Octave [11]. All the computation times given in the paper have been performed on a 2.80 GHz Intel Core i7-4810MQ CPU with 8 GB of memory.

Plan

The structure of this paper is as follows. The class of systems considered and some preliminary definitions are given in Section 2. Our symbolic approach, which is based on the tiling of the state space and backward reachability, is explained in Section 3. In Section 4, we present a centralized method to synthesize a controller based on a "generate-and-test" tiling procedure. A distributed approach is then given in Section 5 where we introduce a state over-approximation technique in order to avoid the use of non-local information by the subsystem controllers. For both methods, we provide attainability and stability guarantees on the controlled trajectories of the system. Finally, in Section 6, our distributed approach is applied on a real-case study of temperature control in a building with 11 rooms and 2 11 = 2048 switching modes of control.

State-dependent Switching Control

Control modes

Consider the discrete-time system with finite control:

x 1 (t + 1) = f 1 (x 1 (t), x 2 (t), u 1 ) x 2 (t + 1) = f 2 (x 1 (t), x 2 (t), u 2 )
where x 1 (resp. x 2 ) is the 1st-component (resp. 2nd-component) of the state vector variable, which takes its values in R n1 (resp. R n2 ), and u 1 (resp. u 2 ) is the 1st (resp. 2nd) component of the control mode variable which takes its values in the finite set U 1 (resp. U 2 ). We will often use x for (x 1 , x 2 ), u for (u 1 , u 2 ), and n for n 1 + n 2 . We will also abbreviate the set U 1 × U 2 as U . Let N be the cardinal of U , and N 1 (resp. N 2 ) the cardinal of U 1 (resp. U 2 ). We have

N = N 1 •N 2 .
More generally, we abbreviate the discrete-time system under the form:

x(t + 1) = f (x(t), u)
where x is a vector state variable which takes its values in

R n = R n1 × R n2 , u is of the form (u 1 , u 2 ) where u 1 takes its values in U 1 and u 2 in U 2 .
In this context, we are interested by the following centralized control synthesis problem: at each discrete-time t, select the appropriate mode u ∈ U in order to satisfy a given property. In this paper we focus on state-dependent control, which means that, at each time t, the selection of the value of u is done by considering only the values of x(t).

In the distributed context, the control synthesis problem consists in selecting concurrently the value of u 1 in U 1 according to the value of x 1 (t) only, and the value of u 2 in U 2 according to the value of x 2 (t) only.

The properties that we consider are attainability properties: given a set S and a set R, we look for a control which will steer any element of S to R in a bounded number of steps. We will also consider stability properties, which means, that once the state x of the system is in R at time t, the control will maintain it in R indefinitely at t + 1, t + 2, • • • . Actually, given a state set R, we will present a method which does not start from a given set S, but constructs it, together with a control which steers all the elements of S to R within a bounded number of steps (S can be seen as a "capture set" of R).

In this paper, we consider that R and S are "rectangles" of the state space. More precisely, R = R 1 × R 2 is a rectangle of reals, i.e., R is a product of n closed intervals of reals, and R 1 (resp. R 2 ) is a product of n 1 (resp. n 2 ) closed intervals of reals. Likewise, we assume that S = S 1 ×S 2 is a rectangular sub-area of the state space.

Example 1. The centralized and distributed approaches will be illustrated by the example of a two rooms apartment, heated by two heaters located in each room (adapted from (Girard)). In this example, the objective is to control the temperature of the two rooms. There is heat exchange between the two rooms and with the environment. The continuous dynamics of the system is given by the equation:

˙ T 1 T 2 = -α 21 -α e1 -α f u 1 α 21 α 12 -α 12 -α e2 -α f u 2 T 1 T 2 + α e1 T e + α f T f u 1 α e2 T e + α f T f u 2 .
Here T 1 and T 2 are the temperatures of the two rooms, and the state of the system corresponds to T = (T 1 , T 2 ). The control mode variable u 1 (respectively u 2 ) can take the values 0 or 1 depending on whether the heater in room 1 (respectively room 2) is switched off or switched on (hence

U 1 = U 2 = {0, 1}). Hence, here n 1 = n 2 = 1, N 1 = N 2 = 2 and n = 2, N = 4.
T e corresponds to the temperature of the environment, and T f to the temperature of the heaters. The values of the different parameters are the following:

α 12 = 5×10 -2 , α 21 = 5×10 -2 , α e1 = 5×10 -3 , α e2 = 5×10 -3 , α f = 8.3×10 -3
, T e = 10 and T f = 35.

We suppose that the heaters can be switched periodically at sampling instants τ, 2τ, . . . (here, τ = 5s). By integration of the continuous dynamics between t and t + τ , the system can be easily put under the desired discrete-time form:

T 1 (t + 1) = f 1 (T 1 (t), T 2 (t), u 1 ) T 2 (t + 1) = f 2 (T 1 (t), T 2 (t), u 2 ),
where f 1 and f 2 are affine functions.

Given an objective rectangle for T = (T 1 , T 2 ) of the form R = [18.5, 22] × [18. [START_REF] Girard | Low-complexity switching controllers for safety using symbolic models[END_REF]22], the control synthesis problem is to find a rectangular capture set S as large as possible, from which one can steer the state T to R ("attainability"), then maintain T within R for ever ("stability").

Control patterns

It is often easier to design a control of the system using several applications of f in a row rather than using just a single application of f at each time. We are thus led to the notion of "macro-step", and "control pattern". A (control) pattern π = (π 1 , π 2 ) of length k is a sequence of modes defined recursively by:

1. π is of the form (u 1 , u 2 ) ∈ U 1 × U 2 if k = 1, 2. π is of the form (u 1 • π 1 , u 2 • π 2 ), where u 1 (resp. u 2 ) is in U 1 (resp. U 2 ), and (π 1 , π 2 ) is a (control) pattern of length k -1 if k ≥ 2.
The set of patterns of length k is denoted by

Π k (for length k = 1, Π 1 = U ). Likewise, for k ≥ 1, we denote by Π k 1 (resp. Π k 2 ) the set of sequences of k ele- ments of U 1 (resp. U 2 ).
For a system defined by

x(t+1) = f (x(t), (u 1 , u 2 )) and a pattern π = (π 1 , π 2 ) of length k, one can define recursively x(t + k) = f (x(t), (π 1 , π 2 )) with (π 1 , π 2 ) ∈ Π k , by: 1. f (x(t), (π 1 , π 2 )) = f (x(t), (u 1 , u 2 )), if (π 1 , π 2 ) is a pattern of length k = 1 of the form (u 1 , u 2 ) ∈ U , 2. f (x(t), (π 1 , π 2 )) = f (f (x(t), (π 1 , π 2 )), (u 1 , u 2 )), if (π 1 , π 2 ) is a pattern of length k ≥ 2 of the form (u 1 • π 1 , u 2 • π 2 ) with (u 1 , u 2 ) ∈ U and (π 1 , π 2 ) ∈ Π k-1 .
One defines (f (x, π)) 1 ∈ R n1 and (f (x, π)) 2 ∈ R n2 to be the 1st and 2nd com-

ponents of f (x, π) ∈ R n1 × R n2 = R n , i.e: f (x, π) = ((f (x, π)) 1 , f (x, π) 2 ).
In the following, we suppose that K ∈ N is an upper bound of the length of patterns. The value of K can be seen as a maximum number of time steps, for which we compute the future behavior of the system ("horizon"). We denote by

Π ≤K 1 (resp. Π ≤K 2 ) the expression 1≤k≤K Π k 1 (resp. 1≤k≤K Π k 2 )
. Likewise, we denote by Π ≤K the expression 1≤k≤K Π k .

3 Control Synthesis Using Tiling

3.1 Tiling Let R = R 1 × R 2 be a rectangle. We say that R is a (finite rectangular) tiling of R if R is of the form {r i1,i2 } i1∈I1,i2∈I2
, where I 1 and I 2 are given finite sets of positive integers, each r i1,i2 is a sub-rectangle of R of the form r i1 × r i2 , and r i1 , r i2 are closed sub-intervals of R 1 and R 2 respectively. Besides, we have

i1∈I1 r i1 = R 1 and i2∈I2 r i2 = R 2 (Hence R = i1∈I1,i2∈I2 r i1,i2
). We will refer to r i1 , r i2 and r i1,i2 as "tiles" of R 1 , R 2 and R respectively. The same notions hold for rectangle S.

In the centralized context, given a rectangle R, the macro-step (backward reachability) control synthesis problem with horizon K consists in finding a rectangle S and a tiling S = {s i1,i2 } i1∈I1,i2∈I2 of S such that, for each (i 1 , i 2 ) ∈ I 1 × I 2 , there exists π ∈ Π ≤K such that:

f (s i1,i2 , π) ⊆ R (i.e., for all x ∈ s i1,i2 : f (x, π) ∈ R).
This is illustrated in Figure 1.

Parametric extension of tiling

In the following, we assume that the set S we are looking for is a parametric extension of R, denoted by R + (a, a), which is defined in the following.

Suppose that We now consider that S is a (parametric) superset of R of the form R+(a, a). We define a tiling S = S 1 × S 2 of S of the form {s i1 × s i2 } i1∈I1,i2∈I2 , which is

R = R 1 × R 2 is given as well as a tiling R = R 1 × R 2 = {r i1 ×r i2 } i1∈I1,i2∈I2 = {r i1,i2 } i1∈I1,i2∈I2 .
obtained from R = R 1 × R 2 = {r i1 × r i2 } i1∈I1,i2∈I2

by a simple extension, as follows:

A tile r i1 (resp. r i2 ) of R 1 (resp. R 2 ) in "contact" with ∂R 1 (resp. ∂R 2 ) is prolonged as a tile s i1 (resp. s i2 ) in order to be in contact with ∂(R 1 + a) (resp. ∂(R 2 + a)); a tile "interior" to R 1 (i.e., with no contact with ∂R 1 ) is kept unchanged, and coincides with s i1 , and similarly for R 2 .

We denote the resulting tiling S by R + (a, a). We also denote s i1 (resp. s i2 ) as r i1 + a (resp. r i2 + a) even if r i1 (resp. r i2 ) is "interior" to R 1 (resp. R 2 ). Likewise, we will denote s i,j as r i,j + (a, a). Note that a tiling of R of index set I 1 × I 2 induces a tiling of R + (a, a) with the same index set I 1 × I 2 , hence the same number of tiles as R, for any a ≥ 0. This is illustrated in Figure 2, where the tiling of R is represented with black continuous lines, and the extended tiling of R + (a, a) with red dashed lines.

Generate and test tilings

By replacing S with R + (a, a) in the notions defined in Section 3.1 the problem of macro-step control synthesis can now be reformulated as:

Finding a tiling R of R which induces a macro-step control of R + (a, a) towards R, for some a ≥ 0; besides, if we find such R, we want to compute the maximum value of a for which the induced control exists. This problem can be solved by a simple "generate and test" procedure: one generates a candidate tiling, then one tests if it satisfies the control property (the control test procedure is explained in Section 4.1); if the test fails, one generate another candidate, and so on iteratively.

In practice, the generation of a candidate R is done, starting from the trivial tiling (made of one tile equal to R), then using successive bisections of R until, either the control test succeeds ("success"), or the depth of bisection of the new candidate is greater than a given upper bound D ("failure"). See details in Appendix 8.1.

Remark 1. Note that, if the generate-and-test process stops with "success" for a tiling R, then the tiling R D,unif orm also solves the problem, where R D,unif orm is the "finest" tiling obtained by bisecting D times all the n components of R. Since R D,unif orm has exactly 2 nD tiles, it is in general impractical to perform directly the control test on it. From a theoretical point of view however, it is convenient to suppose that R = R D,unif orm for reducing the worst case time complexity of the control synthesis procedure to the complexity of the control test part only (see Section 4.1).

Centralized control

Tiling test procedure

As seen in Section 3.2, the (macro-step) control synthesis problem with horizon K consists in finding (the maximum value of) a ≥ 0, and a tiling

R = {r i1,i2 } i1∈I1,i2∈I2 of R such that, for each (i 1 , i 2 ) ∈ I 1 × I 2 , there exists some π ∈ Π ≤K with f (r i1,i2 + (a, a), π) ⊆ R.
In order to test if a tiling candidate R = {r i1,i2 } i1∈I1,i2∈I2 of R satisfies the desired property, we define, for each (i 1 , i 2 ) ∈ I 1 × I 2 :

Π ≤K i1,i2 = {π ∈ Π ≤K | f (r i1,i2 , π) ⊆ R}.
When Π ≤K i1,i2 = ∅, we define:

a i1,i2 = max π∈Π ≤K i 1 ,i 2 max{a ≥ 0 | f (r i1,i2 + (a, a), π) ⊆ R} π i1,i2 = argmax π∈Π ≤K i 1 ,i 2 max{a ≥ 0 | f (r i1,i2 + (a, a), π) ⊆ R} A = min (i1,i2)∈I1×I2 {a i1,i2 }.
For each tile r i1,i2 of R and each π ∈ Π ≤K , the test of inclusion f (r i1,i2 , π) ⊆ R can be done in time polynomial in n when f is affine. Hence the test

Π ≤K i1,i2 = ∅ can be done in O(N K • n α ) since Π ≤K contains O(N K ) elements.
The computation of max{a ≥ 0 |f (r i1,i2 + (a, a), π) ⊆ R} can be done by linear programming in time polynomial in n, the dimension of the state space. The computation time of {a i1,i2 } i1∈I,i2∈I2 , π i1,i2 , and A is thus in O(N K •2 nD ), where D is the maximal depth of bisection. Hence the complexity of testing a candidate tiling R is in O(N K • 2 nD ). By Remark 1 above, the complexity of the control synthesis by generate-and-test is also in O(N K • 2 nD ).

We have:

Proposition 1. Suppose that there exists a tiling R = {r i1,i2 } i1∈I1,i2∈I2 of R such that:

∀(i 1 , i 2 ) ∈ I 1 × I 2 Π ≤K i1,i2 = ∅.
Then R induces a macro-step control of horizon K of R+(A, A) towards R with:

∀(i 1 , i 2 ) ∈ I 1 × I 2 : f (r i1,i2 + (A, A), π i1,i2 ) ⊆ R
where A and π i1,i2 are defined as above.

Once a candidate tiling R satisfying the control test property is found, the generate-and-test procedure ends with success (see Section 3.3), and a set S = R + (a (1) , a (1) ) with a (1) = A has been found. One can then iterate the "generate and test" procedure in order to construct an increasing sequence of nested rectangles of the form R + (a (1) , a (1) ), R + (a (1) + a (2) , a (1) + a (2) ), . . . , which can all be driven to R, as explained in Appendix 8.2.

Example 2. Consider the specification of a two-rooms appartment given in Example 1. Set R = [18.5, 22] × [18. [START_REF] Girard | Low-complexity switching controllers for safety using symbolic models[END_REF]22]. Let D = 1 (the depth of bisection is at most 1), and K = 4 (the maximum length of patterns is 4). We look for a centralized controller which will steer the rectangle S = [18.5-a, 22]×[18.5-a, 22] to R with a as large as possible, and stay in R indefinitely. Using our implementation, the computation of the control synthesis takes 4.14s of CPU time.

The method iterates successfully 15 times the macro-step control synthesis procedure. We find S = R + (a, a) with a = 53.5, i.e. S = [-35, 22] × [-35, 22]. This means that any element of S can be driven to R within 15 macro-steps of length (at most) 4, i.e., within 15 × 4 = 60 units of time. Since each unit of time is of duration τ = 5s, any trajectory starting from S reaches R within 60 × 5 = 300s. Once the trajectory x(t) is in R, it returns in R every macro-step of length (at most) 4, i.e., every 4 × 5 = 20s.

These results are consistent with the simulation given in Figure 3 for the time evolution of (T 1 , T 2 ) starting from (12, 12). Simulations of the control, starting from (T 1 , T 2 ) = (12, 12), (T 1 , T 2 ) = (12, 19) and (T 1 , T 2 ) = (22, 12) are also given in the state space plane in Figure 3.

Stability as a special case of reachability

Instead of looking for a set of the form S = R + (a, a) from which R is reachable via a macro-step, let us consider the particular case where S = R (i.e., a = 0).

The problem is now to construct a tiling R = {r i1,i2 } i1∈I1,i2∈I2 of R such that, for all (i 1 , i 2 ) ∈ I 1 ×I 2 , there exists a pattern π i1,i2 ∈ Π ≤K with f (r i1,i2 , π i1,i2 ) ⊆ R. If such a tiling R exists, then x(t) ∈ R implies x(t + k) ∈ R for some k ≤ K. 6Actually, we can slightly modify the procedure in order to impose, additionally, that ∀k ≤ K x(t + k) ∈ R + ε for some ε > 0 (see Section 5.2). It follows that R + (ε, ε) is stable under the control induced by R.

We can thus treat the stability control of R as a special case of reachability control. 5 Distributed control

Background

In the distributed context, given a set R = R 1 × R 2 , the (macro-step) distributed control synthesis problem with horizon K consists in finding (the maximum value of) a ≥ 0, and a tiling

R 1 = {r i1 } i1∈I1 of R 1 which induces a (macro-step) control on R 1 + a, a tiling R 2 = {r i2 } i2∈I2 which induces a (macro-step) control on R 2 + a.
More precisely, we seek tilings R 1 and R 2 such that: there exists ∈ N such that, for each i 1 ∈ I 1 there exists a sequence π 1 of modes in U 1 , and for each i 2 ∈ I 2 , a sequence π 2 of modes in U 2 such that:

(f ((r i1 +a)×(R 2 +a), (π 1 , π 2 ))) 1 ⊆ R 1 ∧ (f ((R 1 +a)×(r i2 +a), (π 1 , π 2 ))) 2 ⊆ R 2 .
In order to synthesize a distributed strategy where the control pattern π 1 is determined only by i 1 (regardless of the value of i 2 ), and the control pattern π 2 only by i 2 (regardless of the value of i 1 ), we now define an over-approximation X i1 (a, π 1 ) for (f ((r i1 + a) × (R 2 + a), (π 1 , π 2 ))) 1 , and an over-approximation X i2 (a, π 2 ) for (f ((R 1 + a) × (r i2 + a), (π 1 , π 2 ))) 2 . The correctness of these overapproximations relies on the existence of a fixed positive value for parameter ε. Intuitively, ε represents the width of the additional margin (around R + (a, a)) within which all the intermediate states lie when a macro-step is applied to a point of R + (a, a).

Tiling test procedure

Let π k 1 (resp.π k 2 ) denote the prefix of length k of π 1 (resp.π 2 ), and π 1 (k) (resp. π 2 (k)) the k-th element of sequence π 1 (resp. π 2 ). Definition 1. Consider an element r i1 (resp. r i2 ) of a tiling R 1 (resp. R 2 ) of R 1 (resp. R 2 ), and a sequence π

1 ∈ Π ≤K 1 (resp. π 2 ∈ Π ≤K 2 ) of length 1 (resp. 2 ). The approximate 1st-component (resp. 2nd-component) sequence {X k i1 (a, π 1 )} 0≤k≤ 1 (resp. {X k i2 (a, π 2 )} 0≤k≤ 2 ) is defined as follows: -X 0 i1 (a, π 1 ) = r i1 + a and -X k i1 (a, π 1 ) = f 1 (X k-1 i1 (a, π 1 ), R 2 + a + ε, π 1 (k)) for 1 ≤ k ≤ 1 (resp. -X 0 i2 (a, π 2 ) = r i2 + a and -X k i2 (a, π 2 ) = f 2 (R 1 + a + ε, X k-1 i2 (a, π 2 ), π 2 (k)) for 1 ≤ k ≤ 2 ).
We define the property P rop(a, i 1 , π 1 ) of {X k i1 (a, π 1 )} 0≤k≤ 1 by:

-X k i1 (a, π 1 ) ⊆ R 1 + a + ε for 1 ≤ k ≤ 1 -1, and -X 1 i1 (a, π 1 ) ⊆ R 1 .
Likewise, we define the property P rop(a, i 2 , π 2 ) of {X k i2 (a, π 2 )} 0≤k≤ 2 by:

-X k i2 (a, π 2 ) ⊆ R 2 + a + ε for 1 ≤ k ≤ 2 -1, and -X 2 i2 (a, π 2 ) ⊆ R 2 . Figure 4 illustrates property P rop(a, i 1 , π 1 ) for π 1 = (u 1 • v 1 • w 1 ), 1 = 3 and i 1 = 1: in the upper part, P rop(a, i 1 , π 1 ) is not satisfied because X 1 1 (a, π 1 ) ⊆ R 1 + a + ε is false (X 2
1 (a, π 1 ) ⊆ R 1 + a + ε and X 3 1 (a, π 1 ) ⊆ R 1 are true); in the lower part, P rop(a, i 1 , π 1 ) is satisfied. Given a tiling R 1 = {r i1 } i1∈I1 of R 1 , we now define, for each i 1 ∈ I 1 , and each k ∈ {1, . . . , K}:

Π k i1 = {π 1 ∈ Π k 1 | P rop(0, i 1 , π 1 )}.
When Π k i1 = ∅, we define:

a k i1 = max π1∈Π k i 1 max{a ≥ 0 | P rop(a, i 1 , π 1 )} π k i1 = argmax π1∈Π k i 1 max{a ≥ 0 | P rop(a, i 1 , π 1 )}
Given R 2 , we define similarly: Π k i2 , a k i2 and π k i2 .

Suppose now, that: (H1) there exists k 1 ∈ {1, . . . , K} such that ∀i 1 ∈ I 1 :

Π k1 i1 = ∅. (H2) there exists k 2 ∈ {1, . . . , K} such that ∀i 2 ∈ I 2 : Π k2 i2 = ∅.
Then we define a k1 1 = min i1∈I1 {a k1 i1 }, a k2 2 = min i2∈I2 {a k2 i2 }, and A = min{a k1 1 , a k2 2 }.

NB1: Given a tiling R = R 1 × R 2 , (H1) means that the points of R 1 + A can be (macro-step) controlled to R 1 using patterns which all have the same length k 1 ; in other terms, all the macro-steps controlling R 1 + A contain the same number k 1 of elementary steps. And symmetrically for (H2). NB2: The determination of an appropriate value for ε is for the moment done by hand, and is the result of a compromise: if ε is too small, then

f 1 (r i1 + a, R 2 + a, u 1 ) ⊆ R 1 + a + ε; if ε is too big, f 1 (X k-1 i1 , R 2 + a + ε, π 1 (k)) ⊆ R 1 + a.
Given a tiling R = R 1 × R 2 of R and a real ε > 0, the problem of existence and computation of k 1 , k 2 , {π k1 i1 } i1∈I1 , {π k2 i2 } i2∈I2 , and A can be solved by linear programming since f 1 and f 2 are affine. Using the same kinds of calculation as in the centralized case (see Section 4.1), one can see that the complexity of testing Π k i1 = ∅ and Π k i2 = ∅ for 1 ≤ k ≤ K, checking (H1)-(H2), generating k 1 , k 2 , A and {π i1 } i1∈I1 , and

{π i2 } i2∈I2 is in O((max(N 1 , N 2 )) K • 2 max(n1,n2)D ).
Hence the complexity of the control test procedure is also in O((max

(N 1 , N 2 )) K • 2 max(n1,n2)D ). Lemma 1. Consider a tiling R = R 1 × R 2 of the form {r i1 × r i2 } (i1,i2)∈I1×I2 .
Let a ≥ 0. We suppose that (H1) and (H2) hold, and that, for all i 1 ∈ I 1 , P rop(a, i 1 , π 1 ) holds for some π 1 ∈ Π k1 1 , and for all i 2 ∈ I 2 , P rop(a, i 2 , π 2 ) holds for some π 2 ∈ Π k2

2 , then we have:

-in case k 1 ≤ k 2 : (f ((r i1 + a, R 2 + a), (π k 1 , π k 2 ))) 1 ⊆ X k i1 (a, π 1 ) ⊆ R 1 + a + ε and (f ((R 1 + a, r i2 + a), (π k 1 , π k 2 ))) 2 ⊆ X k i2 (a, π 2 ) ⊆ R 2 + a + ε, for all 1 ≤ k ≤ k 1 , and (f ((r i1 + a, R 2 + a), (π k1 1 , π k1 2 ))) 1 ⊆ X k1 i1 (a, π 1 ) ⊆ R 1 , -in case k 2 ≤ k 1 : (f ((r i1 + a, R 2 + a), (π k 1 , π k 2 ))) 1 ⊆ X k i1 (a, π 1 ) ⊆ R 1 + a + ε and (f ((R 1 + a, r i2 + a), (π k 1 , π k 2 ))) 2 ⊆ X k i2 (a, π 2 ) ⊆ R 2 + a + ε, for all 1 ≤ k ≤ k 2 , and (f ((R 1 + a, r i2 + a), (π k2 1 , π k2 2 ))) 2 ⊆ X k2 i2 (a, π 2 ) ⊆ R 2 .
The proof of given in Appendix 9. At t = 0, consider a point x(0) = (x 1 (0), x 2 (0)) of R+(A, A), and let us apply concurrently the strategy induced by R 1 on x 1 , and R 2 on x 2 . After k 1 steps, by Lemma 1, we obtain a point x(k

1 ) = (x 1 (k 1 ), x 2 (k 1 )) ∈ R 1 × (R 2 + A + ε).
Then, after k 1 steps, we obtain again a point x(2k 1 ) ∈ R 1 × (R 2 + A + ε), and so on iteratively. Likewise, we obtain points x(k 2 ), x(2k 2 ), . . . which all belong to (R

1 + A + ε) × R 2 . It follows that, after = lcm(k 1 , k 2 ) steps, we obtain a point x( ) which belongs to R 1 × R 2 = R. Theorem 1. Suppose that there is a tiling R 1 = {r i1 } i1∈I1 of R 1 , a tiling R 2 = {r i2 } i2∈I2 of R 2 ,
and a positive real ε such that (H1) and (H2) hold, and let k 1 , k 2 , A be defined as above.

Let = lcm(k 1 , k 2 ) with = α 1 k 1 = α 2 k 2 for some α 1 , α 2 ∈ N.
Then R 1 induces a sequence of α 1 macro-steps on R 1 +A, and R 2 a sequence of α 2 macro-steps on R 2 + A, such that, applied concurrently, we have, for all i 1 ∈ I 1 and i 2 ∈ I 2 :

(f ((r i1 + A) × (R 2 + A), π)) 1 ⊆ R 1 ∧ (f ((R 1 + A) × (r i2 + A), π)) 2 ⊆ R 2 , for some π = (π 1 , π 2 ) ∈ Π where π 1 (resp. π 2 ) is of the form π 1 1 • • • π α1 1 (resp. π 1 2 • • • π α2 2 ) with π i 1 ∈ Π k1 1 for all 1 ≤ i ≤ α 1 (resp. π i 2 ∈ Π k2 2 for all 1 ≤ i ≤ α 2 ).
Besides, for all prefix π of π, we have

(f ((r i1 +A)×(R 2 +A), π )) 1 ⊆ R 1 +A+ε ∧ (f ((R 1 +A)×(r i2 +A), π )) 2 ⊆ R 2 +A+ε.
If (H1)-(H2) hold, there exists a control that steers R+(A, A) to R in steps. Letting R = R + (A, A), it is then possible to iterate the process on R and, in case of success, generate a rectangle R = R + (A , A ) from which R would be reachable in steps, for some A ≥ 0 and ∈ N. And so on, iteratively, one generates an increasing sequence of nested control rectangles, as in Section 4.1.

Example 3. Consider again the specification of a two-rooms appartment given in Example 1. We consider the distributed control synthesis problem where the 1st (resp. 2nd) state component corresponds to the temperature of the 1st (resp. 2nd) room T 1 (resp. T 2 ), and the 1st (resp. 2nd) control mode component corresponds to the heater u 1 (resp. u 2 ) of the the 1st (resp. 2nd) room.

Set [START_REF] Girard | Low-complexity switching controllers for safety using symbolic models[END_REF]22]. Let D = 3 (the depth of bisection is at most 3), and K = 10 (the maximum length of patterns is 10). The parameter ε is set to value 1.5 • C. We look for a distributed controller which steers any temperature state in S = S 1 × S 2 = [18.5 -a, 22] × [18.5 -a, 22] to R with a as large as possible, then maintain it in R indefinitely.

R = R 1 × R 2 = [18.5, 22] × [18.
Using our implementation, the computation of the control synthesis takes 220s of CPU time.

The method iterates 8 times the macro-step control synthesis procedure. We find S = [18.5 -a, 22] × [18.5 -a, 22] with a = 6.5, i. These results are consistent with the simulation given in Figure 5 showing the time evolution of (T 1 , T 2 ) starting from (12, 12). Simulations of the control are also given in the state space plane, in Figure 5, for initial states (T 1 , T 2 ) = (12, 12), (T 1 , T 2 ) = (12, 19) and (T 1 , T 2 ) = (22, 12).

Not surprisingly, the performance guaranteed by the distributed approach (a = 6.5, attainability of R in 400s) are worse than those guaranteed by the centralized approach of Example 2 (a = 53.5, attainability of R in 300s). However, unexpectedly, the CPU computation time in the distributed approach (220s) is here worse than the CPU time of the centralized approach (4.14s). This relative inefficiency is due to the small size of the example.

Case study

This case study, proposed by the Danish company Seluxit, aims at controlling the temperature of an eleven rooms house, heated by geothermal energy.

The continuous dynamics of the system is the following:

d dt T i (t) = n j=1 A d i,j (T j (t) -T i (t)) + B i (T env (t) -T i (t)) + H v i,j .v j (1)
The temperatures of the rooms are the T i . The matrix A d contains the heat transfer coefficients between the rooms, matrix B contains the heat transfer coefficients betweens the rooms and the external temperature, set to T env = 10 • C for the computations. The control matrix H v contains the effects of the control on the room temperatures, and the control variable is here denoted by v j . We have v j = 1 (resp. v j = 0) if the heater in room j is turned on (resp. turned off). We thus have n = 11 and N = 2 11 = 2048 switching modes.

Note that the matrix A d is parametrized by the open of closed state of the doors in the house. In our case, the average between closed and open matrices was taken for the computations. The exact values of the coefficients are given in [START_REF] Kim G Larsen | Online and Compositional Learning of Controllers with Application to Floor Heating[END_REF]. The controller has to select which heater to turn on in the eleven rooms. Due to a limitation of the capacity supplied by the geothermal device, the 11 heaters cannot be turned on at the same time. In our case, we set to 4 the maximum number of heaters turned on at the same time.

We consider the distributed control synthesis problem where the 1st (resp. 2nd) state component corresponds to the temperatures of rooms 1 to 5 (resp. 6 to 11), and the 1st (resp. 2nd) control mode component corresponds to the heaters of rooms 1 to 5 (resp. 6 to 11). Hence n 1 = 5, n 2 = 6, N 1 = 2 5 , N 2 = 2 6 . We impose that at most 2 heaters are switched on at the same time in the 1 st sub-system, and at most 2 in the 2 nd sub-system.

Let D = 1 (the depth of bisection is at most 1), and K = 4 (the maximum length of patterns is 4). The parameter ε is set to value 0.5 • C. The sampling time is τ = 15 min.

We look for a distributed controller which steers any temperature state in the rectangle S = [18 -a, 22] 11 to R = [18, 22] 11 with a as large as possible, then maintain the temperatures in R indefinitely.

Using our implementation, the computation of the control synthesis takes around 20 hours of CPU time.

The method iterates 15 times successfully the macro-step control synthesis procedure. We find S = [18 -a, 22] 11 with a = 4.2, i.e. S = [13. [START_REF] Liberzon | Switching in systems and control[END_REF]22] 11 . This means that any element of S can be driven into R within 15 macro-steps of length (at most) 4, i.e., within 15 × 4 = 60 units of time. Since each unit of time is of duration τ = 15 min, any trajectory starting from S attains R within 60 × 15 = 900 min. The trajectory is then guaranteed to stay in R + (ε, ε) = [17.5, 22.5] 11 .

These results are consistent with the simulation given in Figure 6 showing the time evolution of the temperature of the rooms, starting from 14 11 .

Robustness simulations for our controller are given in Appendix 10.

Final Remarks

In this paper, we have proposed a distributed approach for control synthesis and applied it to a real floor heating system. To our knowledge, this is the first time that attainability and stability properties are guaranteed for a case study of this size. The method can be extended to take into account obstacles and Let us now explain how we find a tiling R of R such that Π i1,i2 = ∅. We focus on the centralized case, but the distributed case is similar. We start from the trivial tiling R 0 = {R} which consists of just one tile equal to R. If f (R, π) ⊆ R for some π ∈ Π ≤K , then R 0 is the desired tiling. Otherwise, we refine R 0 by bisection, which gives a tiling R 1 of the form {r (i,1),(j,2) } 1≤i,j≤n . If, for all 1 ≤ i, j ≤ n there exists some π ∈ Π ≤K such that f (r (i,1),(j,2) , u) ⊆ R, then R 1 is the desired tiling. Otherwise, there exist some "bad" tiles of the form r (i,1),(j,2) with 1 ≤ i, j ≤ n such that ∀π ∈ Π ≤K f (r (i,1),(j,2) , π) ⊆ R; we then transform R 1 into R 2 by bisecting all the bad tiles. And so on iteratively, we produce tilings R 1 , R 2 , • • • , R d until either no bad tiles remain in R d (success) or d is greater than the upper bound D of depth of bisection (failure).

Iterated macro-step control synthesis

Suppose that we are given an objective rectangle R = R 1 × R 2 . If the one-step control synthesis described in Section 8.1 succeeds, then there is a positive null real a (1) = A and a tiling R of R which induces a control steering all the points of R (1) = R + (a (1) , a (1) ) to R in one step. Now the macro-step control synthesis can be reapplied to R (1) . If it succeeds again, then it produces a tiling R (1) of R (1) which induces a control that steers R (2) = R (1) + (a (2) , a (2) ) to R (1) for some a (2) ≥ 0. And so on, the iterated application of macro-step control synthesis outputs a sequence of tilings R (i) which induce a control that steers R (i+1) = R + (Σ i+1 j=1 a (j) , Σ i+1 j=1 a (j) ) to R (i) , for some a (j) ≥ 0 (1 ≤ j ≤ i + 1). We thus synthesize a control which steers R (i+1) to R in at most i + 1 macro-steps (i ≥ 0), using an increasing sequence of nested rectangles around R. This is illustrated in Figure 7, for i = 1.

The iteration process halts at step, say m, when the last macro-step control synthesis fails because the maximum bisection depth D is reached while "bad" tiles still remain (see Section 8.1). We also stop the process when the last macrostep control synthesis outputs a real a (m) which is smaller than a given bound η > 0: this is because the sequence of controllable rectangles around R seems to approach a limit. Fig. 7. Iterated control of R (1) = R + (a (1) , a (1) ) towards R, and R (2) = R (1) + (a (2) , a (2) ) towards R (1) .

Note that the first argument of f 1 in the last expression satisfies [f ((r i1 + a, R 2 + a), (π k-1 1 , π k-1

2

))] 1 ⊆ X k i1 by P 1 i1 (k -1). Besides, the second argument satisfies [f ((r i1 +a, R 2 +a), (π k-1

1 , π k-1 2 ))] 2 ⊆ j2∈I2 X k-1 j2 ⊆ R 2 +a+ε, because 1. r i1 + a ⊆ R 1 + a 2. j2∈I2 X k-1 j2 ⊆ R 2 + a + ε since X k-1 j2
⊆ R 2 + a + ε (by P 2 j2 (k -1) which holds for all j 2 ) 3. [f ((R 1 + a, r j2 + a), (π k-1 1 , π k-1

2

))] 2 ⊆ X k-1 j2 (by P 1 j2 (k -1)).

Hence

f 1 ([f ((r i1 + a, R 2 + a), (π k-1 1 , π k-1 2
))] 1 , [f ((r i1 + a, R 2 + a), (π k-1

1 , π k-1 2
))] 2 , π

(k) 1 ) ⊆ f 1 (X k-1 i1 , R 2 + a + ε, π 1 (k)) = X k i1
We have thus proved P 1 i1 (k):

(f (r i1 + a, R 2 + a, (π k 1 , π k 2 ))) 1 ⊆ X k i1
This completes the proof of ∀i 1 ∈ I 1 , P 1 i1 (k) ∧ P 2 i1 (k) We prove ∀i 2 ∈ I 2 , P 1 i2 (k) ∧ P 2 i2 (k) similarly, which achieves the proof of P (k). The proof of (f ((r i1 + a, R 2 + a), (π k1 1 , π k1 2 ))) 1 ⊆ X k1 i1 (a, π 1 ) ⊆ R 1 is similar.

10 Appendix: Robustness Experiments

We now perform the same simulations as in Figure 6, except that the environment temperature is not fixed at 10 • C but follows scenarios of soft winter (Figure 8) and spring (Figure 9). The environment temperature is plotted in green in the figures. The spring scenario is taken from [START_REF] Kim G Larsen | Online and Compositional Learning of Controllers with Application to Floor Heating[END_REF], and the soft winter scenario is the winter scenario of [START_REF] Kim G Larsen | Online and Compositional Learning of Controllers with Application to Floor Heating[END_REF] with 5 additional degrees. We see that our controller, which is designed for T env = 10 • C still satisfies the properties of attainability and stability. These simulations are very close those obtained in [START_REF] Kim G Larsen | Online and Compositional Learning of Controllers with Application to Floor Heating[END_REF].

Fig. 1 .

 1 Fig. 1. Mapping of tile s2,3 to R via pattern π2,3, and mapping of tile s3,1 via π3,1.

  R 1 can be seen as a product of n 1 closed intervals of the form [ , m]. Consider a non negative real parameter a. Let (R 1 +a) denote the corresponding product of n 1 intervals of the form [ -a, m + a]. 5 We define (R 2 + a) similarly. Finally, we define R + (a, a) as (R 1 + a) × (R 2 + a).

Fig. 2 .

 2 Fig. 2. Tiling of R + (a, a) induced by tiling R of R.

Fig. 3 .

 3 Fig. 3. Simulations of the centralized reachability controller for three different initial conditions plotted in the state space plane (left); simulation of the centralized reachability controller for the initial condition (12, 12) plotted within time (right).

Fig. 4 .

 4 Fig. 4. Illustration of ¬P rop(a, i1, π1) with i1 = 1, |π1| = 1 = 3 (left), and P rop(a, i1, π1) (rigth).

  e. S = [12, 22] × [12, 22]. This means that any element of S can be driven to R within 8 macro-steps of length (at most) 10, i.e., within 8 × 10 = 80 units of time. Since each unit of time is of duration τ = 5s, any trajectory starting from S reaches R within 80 × 5 = 400s. The trajectory is then guaranteed to always stay (at each discrete time t) in R + (ε, ε) = [17, 23.5] × [17, 23.5].

Fig. 5 .

 5 Fig. 5. Simulations of the distributed reachability controller for three different initial conditions plotted in the state space plane (left); simulation of the distributed reachability controller for the initial condition (12, 12) plotted within time (right).

Fig. 8 .

 8 Fig. 8. Simulation of the Seluxit case study in the soft winter scenario.

Fig. 9 .

 9 Fig. 9. Simulation of the Seluxit case study in the spring scenario.

This separability technique is made possible because the difference equation x1(t + 1) = f1(x1(t), x2(t), u1) (see Section 2.1) does not involve the control mode u2.

Actually, we will consider in the examples that (R1 + a) is a product of intervals of the form [ -a, m] where the interval is extended only at its lower end, but the method is strictly identical.

If x(t) ∈ R, then x(t) ∈ ri,j for some (i, j) ∈ I1 × I2, hence x(t + k) = f (x, πi,j) ∈ R for some k ≤ K.

safety constraints. We are currently investigating an extension of the method to systems with non linear dynamics and varying parameters. 9 Appendix: Proof of Lemma 1

and by P 2 i1 (k) X k i1 ⊆ R 1 + a + ε and similarly for P 1 i2 (k) and P 2 i2 (k). We are going to show by induction on k the property P (k):

1 and π k 2 are of the form u 1 and u 2 . We have:

Let us now suppose that k > 1 and that P (k -1) holds. Let us prove P (k). Properties P 2 i1 (k) and P 2 i2 (k) are true for all i 1 , i 2 because, by construction, the sequence X k i1 (resp. X k i2 ) satisfies P rop(a, i 1 , π 1 ) (resp. P rop(a, i 2 , π 2 )). Let us prove P 1 i1 (k) and P 1 i2 (k):

)), (π

))] 1 ,

[f ((r i1 + a, R 2 + a), (π k-1

))] 2 , π 1 (k)).