Z Ayane Chadia
Z Ayane Chadia
Doctorat ParisTech
THÈSE
pour obtenir le grade de docteur délivré par
Chadia ZAYANE
le 11 janvier 2011
Jury
1 Introduction 1
1.1 Motivations et enjeux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.2 Enjeux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Problématique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 Présentation du problème . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.2 Données et contraintes . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.2.1 Présentation des données . . . . . . . . . . . . . . . . . . . . 4
1.2.2.2 Les contraintes . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.3 Objectifs de la thèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Approche proposée . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Organisation du document . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
i
3 Identification et problèmes mal posés 27
3.1 La procédure d’identification . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1.1 Les étapes préliminaires . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1.2 Structures de modèle et étapes clés . . . . . . . . . . . . . . . . . . . . 29
3.1.3 Enjeu et limites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2 Identifiabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2.1 Définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.2.2 Expérience informative pour une structure donnée . . . . . . . . . . . . 35
3.2.2.1 Contenu fréquentiel . . . . . . . . . . . . . . . . . . . . . . 36
3.2.2.2 Erreur de convergence et matrice d’information de Fisher . . 37
3.3 Identification d’un système MISO ou en boucle fermée . . . . . . . . . . . . . 38
3.3.1 Identification d’un système multi entrées . . . . . . . . . . . . . . . . . 38
3.3.2 Identification en boucle fermée . . . . . . . . . . . . . . . . . . . . . . 39
ii
5.2.3.1 Le filtre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.2.3.2 Le lisseur . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.2.3.3 L’échantillonneur . . . . . . . . . . . . . . . . . . . . . . . . 74
Annexes 114
iii
B Matrice d’information de Fisher/ algorithme du gradient 123
B.1 Calcul de la matrice d’information de Fisher . . . . . . . . . . . . . . . . . . . 123
B.1.1 Évaluation de ∇2θ log [ fθ (Y )] pour un système récursif . . . . . . . . . . 123
B.2 Algorithme de gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
Bibliographie 131
iv
Table des figures
v
4.6 Écart entre courbe de charge de référence et simulée, pour une perturbation de
1 1
référence ref_bruit : ΣX2 0 = diag([0.2; 0.1; 0.2]), Q 2 = diag([0.4; 0.3; 0.3]), R = 1 . 50
4.7 Écart quadratique entre courbe de charge de référence et simulée, pour une pertur-
bation dans le sens du gradient, grad : gradient normalisé . . . . . . . . . . . . . 51
4.8 Sensibilité Qch par rapport à une variation de chacun des paramètres . . . . . . . 52
4.9 Données CLIM 2000 pour les 5 zones . . . . . . . . . . . . . . . . . . . . . . . . 53
4.10 Données du bâtiment considéré comme une zone thermique unique . . . . . . . . . 54
4.11 Densité Spectrale de Puissance des différents signaux . . . . . . . . . . . . . . . . 54
4.12 Courbes de charges réelle et générée par le modèle identifié . . . . . . . . . . . . 55
4.13 Distribution de l’erreur de prédiction . . . . . . . . . . . . . . . . . . . . . . . . 56
4.14 Autocorrélation des résidus et leur intercorrélation avec Text . . . . . . . . . . . . 57
4.15 Intercorrélation entre les résidus et les entrées Qs et E . . . . . . . . . . . . . . . 58
vi
6.14 Simulations de l’état associé aux données CLIM2000, chaque couleur correspond à
une itération de l’algorithme de Gibbs, la température intérieure moyenne fournie
par CLIM2000 est reportée sur le premier graphe (noir). . . . . . . . . . . . . . . 94
6.15 Histogrammes des distributions a priori (violet) et a posteriori (vert) du gain sta-
tique par rapport au flux solaire . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
6.16 Données recueillies au bâtiment LGEP1 . . . . . . . . . . . . . . . . . . . . . . . 96
6.17 Histogrammes des distributions a posteriori des paramètres du modèle du bâtiment
LGEP1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
s du bâtiment LGEP1, température moyenne
6.18 Simulations de la température intérieure Tint
Tint enregistrée (noir). L’écart entre chaque simulation et la moyenne est représenté
sur le deuxième graphe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6.19 Température intérieure d’un bureau appartenant au bâtiment LGEP1, pour une pé-
riode semblable (hiver) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
6.20 Saturation de la commande chauffage . . . . . . . . . . . . . . . . . . . . . . . . 100
6.21 Effet de la saturation sur le comportement du système en boucle fermée . . . . . . 101
6.22 Schéma bloc de la boucle fermée avec anti windup . . . . . . . . . . . . . . . . . 101
6.23 Températures et puissance de chauffage en fonction du gain Ks . . . . . . . . . . . 102
6.24 Régime transitoire avec/sans saturation . . . . . . . . . . . . . . . . . . . . . . . 103
6.25 Comparaison du taux d’acceptation en fonction du pas de discrétisation . . . . . . 113
6.26 Exemple de trajectoires hamiltoniennes . . . . . . . . . . . . . . . . . . . . . . . 113
vii
Liste des tableaux
viii
Liste des algorithmes
ix
Liste des symboles
xi
Chapitre 1
Introduction
1.1.1 Motivations
En France, le secteur du bâtiment (résidentiels, tertiaires et industriels) est le secteur ayant
la plus grande part dans la demande énergétique parmi tous les secteurs économiques. En ef-
fet, ses besoins représentent environ 43% de l’énergie finale 1 totale, ce qui correspond à une
consommation annuelle autour de 400 kWh d’énergie primaire par m2 chauffé 2 . Cette part de
consommation du bâtiment ne cesse d’augmenter, par exemple les logements et les bureaux ont
vu leur consommation augmenter de +30% durant les 30 dernières années. Cette augmentation
est liée :
ä à la croissance permanente de la taille du parc considéré ;
ä à l’augmentation de la surface moyenne occupée ;
ä à l’accroissement du confort ;
ä et à l’émergence de plusieurs nouveaux usages électriques.
Dans le cas particulier du secteur résidentiel, le chauffage électrique représente la plus grande
part de consommation (supérieure à 60% de la consommation totale). Une meilleure gestion
de la consommation électrique du chauffage devient une urgence, notamment dans un contexte
économique qui évolue rapidement et face à une conscience accrue de la contrainte environ-
nementale. Concrètement, il s’agirait de choisir une installation adaptée aux besoins (gestion
de la demande en énergie) et de s’assurer de la bonne isolation du bâtiment, qu’il s’agisse de
bâtiments existants ou en cours de construction.
1. L’énergie finale est l’énergie distribuée à tous les consommateurs
2. Données ADEME 2005
1
1.1.2 Enjeux
Une bonne connaissance du comportement thermique du bâtiment (résidentiel, tertiaire, indus-
triel) contribue à améliorer la gestion de la demande en énergie du chauffage électrique. Le
commercialisateur exploite cette information pour proposer des offres de services énergétiques
tels que le diagnostic du bâtiment et des recommandations d’amélioration de la structure du
bâtiment ou de remplacement du système de chauffage électrique. Pour le client, une meilleure
connaissance du comportement de son bâtiment est nécessaire pour modifier son comportement
énergétique, pour réduire sa facture ou encore pour améliorer son confort.
1.2 Problématique
1.2.1 Présentation du problème
Le comportement thermique d’un bâtiment dépend de différents facteurs illustrés sur la figure
1.1. Ces derniers peuvent être classés en trois familles. La première famille englobe les facteurs
intrinsèques à la structure du bâtiment lui-même (composition, surface à chauffer, etc.). La
seconde famille regroupe les données météorologiques 3 ayant une influence sur le bâtiment
tels que la température extérieure et le flux solaire. La troisième famille, dont l’interaction avec
le bâtiment est expliquée dans le paragraphe suivant, regroupe les facteurs « contrôlables » par
le client, à savoir la puissance appelée par le système de chauffage électrique, la température de
consigne choisie et enfin les apports internes (apports calorifiques autre que le chauffage) dont
les apports gratuits dus à la présence de personnes.
Apports
internes
Flux Température
solaire de consigne
Bâtiment
Température Consommation
extérieure chauffage
Température
intérieure
Le système de chauffage a pour rôle d’ajuster la température intérieure afin de maintenir le ni-
veau de confort souhaité. De ce fait, il constitue la partie « pilotable » du bâtiment. Cet usage
électrique fournit la puissance électrique nécessaire à la réduction de l’écart entre la tempéra-
ture de consigne et la température ambiante. Ainsi, la consommation du chauffage devient une
3. Nous utiliserons désormais le mot « météo » pour désigner les données météorologiques
2
commande pour le bâtiment. Il faudrait donc prendre en compte cette puissance, qui renseigne
sur le système de régulation et sur l’état du bâtiment, dans le modèle global. Ce couplage entre
le bâtiment et le système de chauffage traduit le lien entre le confort souhaité par les usagers et
la réponse thermique du bâtiment. Soulignons que l’action du système de chauffage transforme
le fonctionnement du bâtiment d’un système en boucle ouverte (bâtiment seul) en un système en
boucle fermée (la sortie agit sur l’entrée par rétroaction).
Les apports internes peuvent être vus comme une « perturbation » du comportement thermique
du bâtiment, étant donné que cet apport calorifique vient s’ajouter à l’apport calorifique généré
par le chauffage électrique injecté dans le bâtiment.
On aboutit ainsi à une première modélisation en bloc du comportement du bâtiment prenant en
compte les éléments décrits précédemment présentée sur la figure 1.2.
Courbe de charge
de chauffage
Température
extérieure
Flux
Apports solaire
internes
Température Température
de consigne +- Régulation ++ Bâtiment intérieure
F IGURE 1.2: Schéma bloc associé au modèle global du bâtiment : entrées (vert), mesures (noir)
et systèmes inconnus (rouge)
Sur cette figure, les entrées du système de chauffage sont la température de consigne, grandeur
choisie par l’usager et souvent inconnue et la température intérieure (ambiante). La sortie de
la régulation est une puissance donnée du chauffage, grandeur qui peut être mesurée. Cette
puissance électrique est convertie en énergie calorifique injectée dans bâtiment. Par ailleurs,
le fonctionnement de ce dernier est perturbé par les apports internes et les données météo. La
résultante de ces interactions est une température intérieure donnée, qui va à son tour être prise
en compte par le système de régulation du chauffage. La sortie du système est une température
intérieure du bâtiment, grandeur dont la définition s’avère délicate comme nous l’illustrerons
ultérieurement.
Le schéma bloc présenté sur la figure 1.2 met en évidence que les deux composantes du sys-
tème, le chauffage et le bâtiment, sont couplées. La prise en compte de ce couplage a transformé
le problème en régime libre en un problème en boucle fermée. L’identification du modèle glo-
bal, cas qui nous intéresse dans cette thèse, nécessite alors des précautions quant à l’applica-
tion des approches classiques d’identification généralement conçues pour des systèmes opérant
en boucle ouverte. La prise en compte de la régulation augmente également la complexité du
problème, et ce à travers l’introduction de paramètres inconnus (du régulateur) qui viennent
s’ajouter aux inconnus du modèle thermique à proprement parler.
3
1.2.2 Données et contraintes
Dans ce paragraphe, nous décrivons les données d’entrée/sortie du système global « chauffage
+ bâtiment », les connaissances a priori et/ou contraintes sur ces données.
Les données exploitées dans nos travaux proviennent de sources variées et ont des propriétés
assez différentes :
Données électriques La puissance appelée par le chauffage et les apports internes sont obte-
nus par une décomposition de la courbe de charge électrique totale 5 . Disposer de ces grandeurs
implique l’utilisation d’une méthode de décomposition de la CdC (Courbe de Charge).
Un logiciel de décomposition de la CdC (Barrois et Jicquel [5]), appelé « Fdd_Usages », adapté
au secteur tertiaire bureau été développé au sein d’EDF. Il permet de fournir à partir d’une
courbe de charge annuelle échantillonnée au pas de 10 minutes :
1. un profil de journée minimum 6 , ou talon minimum, correspondant à la consommation des
usages permanents correspondant généralement à la journée du dimanche ;
2. un talon type de chaque journée de la semaine correspondant aux usages liés à la présence
de personnes (ordinateurs, éclairage,...) ;
3. une variation saisonnière, composante restante de la courbe de charge, attribuée à la
consommation de chauffage en hiver et de climatisation en été.
Nous soulignons le fait que la part de chauffage obtenue par cette méthode est généralement
légèrement sur-estimée. Cela est du à la présence d’usages dépendant de la saison et non pris
en compte comme l’éclairage extérieur.
4. L’intermittence est un signal sous forme de créneaux, la valeur minimale est appelée « réduit » et la valeur
maximale « confort ». La période des créneaux est généralement journalière.
5. Les cas que nous étudierons concernent des bâtiments où tous les usages sont électriques mais l’approche
est généralisable pour d’autres types d’énergies (gaz, fuel,...).
6. Pour chaque pas de temps, le minimum est pris en considérant les consommations de toutes les journées de
l’année à la même heure.
4
Île-de-France
H1a
Trappes
H1b Nancy
H2a Rennes
H2b
Mâcon
La Rochelle
H1c
Agen H2d
H2c Carpentras
H3 Nice
F IGURE 1.3: Zones climatique définies par la Réglementation Thermique de 2005 (RT2005)
Données météo En pratique, la température extérieure et le flux solaire sont les grandeurs mé-
téo les plus influentes sur le comportement thermique du bâtiment et sont fournies par la station
météo la plus proche. En effet, Météo-France propose une large gamme de données climatolo-
giques avec des pas de temps allant de l’horaire au mensuel. Les observations sont relevées sur
les quelques 550 stations du réseau RADOME de Météo-France en France métropolitaine ainsi
que ceux des partenaires avec un total de 1250 stations, sachant que la France est répartie en
huit zones climatiques comme l’illustre la figure 1.3.
Le découpage définie par la RT2005 permet de distinguer les zones selon le besoin de chauffage
en hiver (H1, H2 et H3) et de confort en été (indice a, b, c et d), la zone H3 étant en zone d’été
d.
Cette réglementation, applicable à partir de septembre 2006, permet de définir en fonction de
la zone climatique, un ensemble de règles à respecter pour toute construction de bâtiment neuf.
Ces règles permettent notamment d’assurer :
ä une consommation d’énergie inférieure à une valeur de référence définie par la RT,
ä un confort d’été en imposant une température intérieure du bâtiment inférieure à une tem-
pérature de référence,
ä le respect de certaines caractéristiques thermiques des éléments du bâtiment.
Les données que nous traiterons sont enregistrées au niveau de la cellule test ETNA 7 du site
des Renardières d’EDF. Les mesures sont faites avec un pas de temps d’une minute et incluent
7. Il s’agit d’un bâtiment conçu à la fin des années 1980 pour tester les performances de différents composants
de bâtiment : composants de paroi, panneaux solaires, systèmes de chauffage/climatisation...
5
F IGURE 1.4: Les différentes composantes du flux solaire, d’après Mazria [46]
plusieurs données comme la vitesse du vent ou la température du ciel. Le flux solaire reçu par
un bâtiment est composé de deux parties comme l’illustre la figure 1.4 :
ä flux diffus (ou diffuse), dépendant de l’état du ciel (présence de nuages,...) ;
ä flux direct : composante dont les apports sont les plus considérables.
La pertinence des données météo dépend de l’emplacement du bâtiment, une maison mitoyenne
par exemple (ou appartement) est moins « sensible » aux influences de la météo qu’une maison
qui n’est pas mitoyenne.
Température intérieure / ambiante Définir une température unique pour un bâtiment donné
pose une difficulté due à l’hétérogénéité des systèmes de chauffage de chaque pièce, à la va-
riabilité des degrés de confort souhaités dans les différentes pièces et à la difficulté de définir
une température d’une pièce notamment en présence de phénomènes irréguliers (ouverture de
fenêtres par exemple).
Considérons l’exemple d’une résidence où chaque pièce est équipée d’un système de chauffage
(convecteur, radiateur, etc.), les systèmes de chauffages élémentaires étant indépendants les uns
des autres et ayant des modes de régulation différents. La chambre et le séjour sont équipés de
sondes de températures. Les profils de température relevés dans la chambre et dans le séjour
présentés sur la figure 1.5 sont très différents et mettent en évidence une différence entre les
températures de confort souhaitées dans les deux pièces (21° dans la chambre et 23° dans la
séjour).
En pratique, c’est le degré de confort de l’occupant qui prime. C’est pourquoi, différentes études
se sont intéressées à la température ressentie. Qu’il s’agisse de température intérieure ou de
température sentie, se pose la question du nombre de zones thermiques à considérer en fonction
de l’uniformité ou non du degré de confort souhaité dans les différentes pièces.
6
Chambre
23
Séjour
Moyenne
22
Température (°C) 21
20
19
18
17
16
15
3900 4000 4100 4200 4300 4400 4500 4600
Temps (10 min)
F IGURE 1.5: Températures intérieures (maison instrumentée par EDF) : a) chambre à coucher,
b) séjour, c) température moyenne
Dans le cadre de la thèse, nous considérons le cas mono-zone : une seule température est définie
en fonction de la température de consigne choisie.
Données « peu riches » pour identifier un modèle du bâtiment Une stratégie d’identifica-
tion d’un modèle donné du bâtiment (sans prendre en compte le système de chauffage) consiste
à exploiter une série d’expériences où l’on choisit les différents scénarios d’excitation adéquats
et les différentes sollicitations. Cette approche est utilisée par exemple pour étudier le com-
portement du système élémentaire composé des parois du bâtiment. La richesse des données
d’entrée permet d’assurer une bonne identification du modèle étudié. Dans le cas étudié dans la
thèse, à savoir l’identification du système global « chauffage + bâtiment », nous n’avons pas la
possibilité de planifier des expériences.
7
D’autre part, selon la gestion du bâtiment, la température de consigne est toujours informative :
ä si elle est irrégulière : le besoin du client n’est pas formulé précisément, il est donc diffi-
cile d’évoquer l’économie d’énergie notamment en minimisant l’énergie consommée par le
chauffage sans se donner une référence fixe.
ä si elle est à scénarios prédéfinis : la plupart des bâtiments, surtout dans le secteur tertiaire,
utilisent des scénarios fixes de gestion du chauffage électrique, il s’agit alors d’exploiter
cette information d’autant plus que la connaissance du bâtiment est limitée.
8
F IGURE 1.6: Modèle électrique équivalent du bâtiment - modèle R3C2
d’identification proposée repose sur le choix d’un modèle thermique du bâtiment qui soit per-
tinent pour l’estimation des constantes temporelles attendues et dont la complexité est réduite
afin que la méthode soit applicable industriellement.
En pratique, le modèle choisi, appelé modèle R3C2 est un modèle électrique équivalent, a été
proposé par EDF et est illustré sur la figure 1.6. C’est un modèle à cinq paramètres, dont les
entrées sont constituées de données météo (température extérieure et flux solaire) et de la courbe
de charge électrique (chauffage et apports internes) et la sortie est la température intérieure du
bâtiment.
Concernant la régulation, nous considérons un régulateur Proportionnel Intégral illustré sur la
figure 1.7. Il s’agit d’un régulateur qui adapte la puissance de chauffe à l’écart entre la tempé-
rature ambiante et la consigne, et qui est largement utilisé dans les systèmes industriels pour sa
simplicité et son efficacité.
Nous soulignons le fait que l’action du régulateur est soumise à une saturation (indépendam-
ment du système de régulation choisi). En effet la puissance de chauffage est positive et infé-
rieure à une valeur maximale Pmax souscrite par l’usager.
L’approche que nous proposons, consiste d’abord à formuler le problème d’identification de
manière à ce que les températures intérieure et de structure (variables intéressantes pour faire le
diagnostic thermique de bâtiment) fassent partie de l’état caché du système et que la courbe de
charge chauffage en soit une observation. Cette formulation aboutit à un système de dynamique
9
linéaire mais qui est non linéaire en les paramètres inconnus.
Par ailleurs, nous avons mené une étude préliminaire d’identification (voir annexe A) avec un
modèle de bâtiment plus simple (d’ordre 1) en se basant sur des connaissances sectorielles et
thermiques fortes. Cette démarche a démontré la nécessité de vérifier la cohérence et la compa-
tibilité entre le modèle choisi, les informations contenues dans les mesures et les connaissances
a priori sur le comportement thermique de bâtiment.
Fort de ce constat, nous avons ensuite choisi une méthode d’estimation, adaptée pour l’identi-
fication du système global et qui s’inscrit dans le cadre classique d’identification paramétrique.
Même si le modèle adopté (R3C2+PI) est « globalement identifiable », c’est à dire que d’un
point de vue structurel, la solution du problème d’identification est unique, la qualité des me-
sures, recueillies en mode de fonctionnement normal du bâtiment sur une période de quelques
jours, ne permet pas de déterminer les paramètres inconnus avec une grande précision.
Ainsi, pour enrichir le contenu informationnel des données, nous avons exploité des connais-
sances a priori sur les paramètres, basées sur une expertise en thermique du bâtiment. La dé-
marche que nous avons proposée consiste à se placer dans le cadre de l’inversion bayésienne en
utilisant l’échantillonneur de Gibbs que l’on conjugue avec le formalisme de lissage de Kalman
pour obtenir des simulations stochastiques des paramètres et des températures. Cette méthode
nous a permis d’aller au delà de la réponse à la problématique par des estimations « ponc-
tuelles », en montrant les causes d’une précision insuffisante et en pointant l’importance du
dosage entre la connaissance a priori sur le modèle de comportement thermique du bâtiment et
les hypothèses sous-jacentes d’un côté et le comportement réel du bâtiment de l’autre.
10
avons consacré une partie de ce chapitre aux approches d’identification en boucle fermée. En
effet, cette dernière requiert la disponibilité de mesures de la grandeur régulée (température
intérieure) et de la grandeur qui assure la régulation (consommation de chauffage).
Le quatrième chapitre est dédié à l’application de la démarche d’identification décrite au cha-
pitre 3 à un jeu de données provenant d’une modélisation fine d’un bâtiment. Ce dernier est
composé de cinq zones thermiques (pièces) et modélisé à l’aide du logiciel CLIM2000 (dé-
veloppé par EDF). Pour identifier son comportement, le bâtiment a été assimilé à une zone
thermique unique avec une consommation électrique, échantillonnée au pas de temps de 5 mi-
nutes et égale à la somme des consommations de toutes les zones. Bien que les résultats de
l’identification soient satisfaisants en ce qui concerne l’aspect prédictif du système identifié, la
précision sur les paramètres suggère l’exploitation d’un élément important : la connaissance a
priori thermique et sectorielle.
Ainsi, le cinquième chapitre nous place dans un cadre plus général de problème inverse, que
nous nous proposons de résoudre sous un angle bayésien. En effet, non seulement les connais-
sances a priori sont plus facilement exploitables mais la méthode, basée sur l’échantillonneur
de Gibbs, que nous proposons fournit des simulations stochastiques des inconnues (incluant
paramètres du modèle global et températures) plus riches en information : elles offrent la pos-
sibilité de construire plusieurs estimateurs, donnent plus de renseignements sur la dispersion de
chaque variable y compris dans sa dissymétrie, avantage qu’un calcul de variance ne permet pas
de réaliser.
Enfin le sixième chapitre est consacré aux résultats de l’inversion bayésienne appliquée au bâ-
timent. L’algorithme est d’abord appliqué à un jeu de données de synthèse afin de tester ses
performances et d’évaluer l’amélioration ou la dégradation des résultats vis à vis de plusieurs
facteurs comme la pertinence de la connaissance a priori sur les paramètres ou la condition
initiale sur les températures. Ce cas test a également servi pour illustrer quelques heuristiques
classiquement utilisées pour diagnostiquer la convergence des algorithmes de simulations sto-
chastiques.
Les données CLIM2000 ont ensuite été étudiées et les résultats obtenus donnent des éléments
de réponse sur la précision obtenue par la méthode d’identification du chapitre 4.
L’application de l’algorithme d’inversion bayésienne sur des données réelles d’un bâtiment
constitué de plusieurs bureaux, en occupation normale, remet en question certaines hypothèses,
notamment l’approximation mono-zone lorsque les consignes de températures dans les différent
bureaux ne sont pas synchrones et la non prise en compte de la saturation. En ce qui concerne
cette difficulté, nous justifions la nécessité de sa prise en compte et proposons quelques pistes
dans le cas où le régulateur est une simple action proportionnelle.
En guise de perspectives à ce travail, nous suggérons la voie de la simulation hamiltonienne
qui, du fait qu’elle ne requiert que des conditions de différentiabilité, pourrait résoudre aussi
bien le problème de la saturation que celui de bâtiments multi-zones (nous pensons en effet à sa
modélisation sous forme de réseau de « zones »).
11
Chapitre 2
Introduction
Ce chapitre présente une synthèse des méthodes de modélisation du comportement thermique,
classées en deux familles : les approches statiques et les approches dynamiques ainsi qu’un
aperçu sur les méthodes d’identification en thermique du bâtiment. Nous détaillons ensuite un
cas particulier de ces modèles, à savoir le modèle R3C2 développé par EDF et qui est le modèle
du bâtiment retenu dans la thèse compte tenu des objectifs et des spécificités du cas d’applica-
tion considéré.
Deux démarches ont été adoptées pour modéliser le bâtiment. La première s’intéressait à une
modélisation simplifiée pour s’affranchir des limitations des moyens de calcul. Les modèles
proposés sont d’ordre réduit et les paramètres sont issus de relevés sur site (puissance électrique,
13
température intérieure, etc.). La deuxième s’intéressait à la compréhension des phénomènes
d’échange thermiques dans le bâtiment à des fins de simulation 1 .
Nous proposons de classer les modèles existants en deux familles en fonction de l’approche
adoptée, à savoir une approche statique et une approche dynamique. Le choix de l’approche est
lié essentiellement à la simplicité du modèle adopté et au pas de temps choisi.
où :
ä Q : puissance du système de chauffage (W ) ;
ä L : coefficient de déperditions statiques globales (W /°C) (dénommé GV en France, UA aux
États Unis) ;
ä Tint : température intérieure moyenne du bâtiment (°C) ;
ä Text : température extérieure moyenne (°C) ;
ä As : surface équivalente sud (m2 ), notée parfois S, est une surface fictive de paroi verticale
orientée au sud recevant la même énergie solaire que l’ensemble des vitrages du bâtiment.
elle est obtenue à partir de la taille des vitrages en prenant en compte un taux de transmission
moyen de 0.6 et en multipliant par un facteur dépendant de l’orientation. Ce facteur vaut
généralement 1 pour le sud, 0.5 pour l’est et l’ouest et 0 pour le nord.
ä I : ensoleillement global vertical sud (W /m2 ), c’est l’énergie solaire reçue par une paroi
verticale orientée au sud ;
ä ε (W ) : grandeur dépendant de l’état des variables mesurées en début et en fin de la période
d’observation, elle est pondérée par l’inverse du pas de mesure 1/∆t.
La différence entre les modèles proposés est liée au choix du nombre et des hypothèses faites
sur les excitations du système étudié, la température extérieure étant un facteur commun à tous
les modèles statiques.
Les modèles existants constituent soit une simplification de l’équation 2.1 en négligeant ε et
éventuellement les apports solaires, soit un enrichissement en développant ε sous forme de
contributions d’autres sollicitations comme la vitesse du vent par exemple. Dans la suite, deux
exemples sont présentés pour illustrer cette famille de modèles.
1. Dans ce contexte, la « simulation » signifie : quantifier la réaction du bâtiment, au sens variation de la
température ou d’énergie consommée, quand toutes les informations sur le bâtiment, les données météo et les
usages sont connues et parfaitement maîtrisées.
2. l’enveloppe d’un bâtiment est l’ensemble des parois qui constituent sa structure.
14
2.1.1.1 Premier exemple : signature énergétique simple
L’équation 2.2 peut encore se simplifier avec l’hypothèse de variations très faibles de la tempé-
rature intérieure du bâtiment, le bilan thermique du bâtiment ne nécessite plus les mesures de la
température intérieure et peut se mettre sous la forme :
Q = α + β Text (2.3)
Les coefficients α et β sont estimés par régression linéaire à partir d’un nuage de point cor-
respondant aux données journalières, hebdomadaires ou mensuelles collectées au niveau du
bâtiment comme l’illustre la figure 2.1.
Q (kW)
7
0 Text (°C)
-5 0 5 10 15 20
La droite de régression obtenue est appelée signature énergétique du bâtiment. Elle a été utili-
sée dès le début des années 80 pour étudier les performances d’une banque de données consti-
tuée de 128 logements résidentiels aux états unis et en Europe (Ribot et al [57]).
3. La même démarche peut être menée pour la climatisation.
15
Variante utilisant le Degré Jour Il peut aussi être intéressant de comparer des bâtiments
dans des zones climatiques différentes ou d’évaluer leurs besoins de chauffage non pas pour
une température extérieure donnée mais en fonction de la sévérité du climat. On définit alors la
notion de Degré Jour (DJ) (Météo France [51]) :
Définition 2.1.1.1. Pour un lieu donné, le Degré Jour est une valeur représentative de l’écart
entre la température d’une journée donnée et un seuil de température préétabli.
Il sert en général à évaluer les dépenses en énergie pour le chauffage ou la climatisation.
Pour un lieu donné, le Degré Jour correspondant au jour J est calculé, à partir des températures
extrêmes de ce même jour, par la méthode dite « Météo » 4 détaillée par l’algorithme 2.1 :
Q (kW)
7
0
0 50 100 150 200 250 300 350
DJ mensuels
4. Il existe une autre méthode de calcul de DJ appelée « méthode Professionnels de l’énergie » plus adaptée
aux sociétés d’exploitation de chauffage
16
2.1.1.2 Deuxième exemple : variation de température extérieure et d’ensoleillement
Ici, l’idée est d’enrichir le modèle précédent et ce en prenant en compte la variation du flux
solaire en plus de l’influence de la température extérieure. La droite de régression linéaire est
décrite par l’équation suivante :
Q = a + bText + cI
Les modèles statiques dits de « signature énergétique » présentés précédemment sont des mo-
dèles simples. Ils offrent la possibilité d’avoir une première caractérisation du bâtiment (pour le
diagnostic par exemple) à travers l’estimation de son gain statique. Mais, ils présentent les trois
limitations suivantes :
ä ils nécessitent une durée d’observation assez longue, allant d’une dizaine de semaines jus-
qu’à plusieurs périodes de chauffe (donc d’années) si l’on travaille avec des pas de temps
mensuels ;
ä ils ne prennent pas en compte le comportement du bâtiment en régime transitoire, compor-
tement qui est utile pour une connaissance plus fine du bâtiment ;
ä ces modèles ne permettent pas de rendre compte de la régulation du chauffage, bien que les
deux composantes interagissent nécessairement.
Ces modèles sont obtenus en enrichissant les modèles statiques de signature énergétique pré-
sentés précédemment afin de prendre en compte la réponse du bâtiment aux différentes solli-
citations au cours du temps. En effet, le modèle statique n’est pas réaliste au sens où il n’est
pas envisageable de substituer les sollicitations par leurs valeurs moyennes. Pour illustrer cette
17
famille de modèles, nous citons le modèle développé par EDF (Gharbi et al [25]) et décrit par
l’équation suivante :
18
Conduction Convection Rayonnement
ä Rayonnement : pour ce mode de transfert, il n’y a pas de contact entre le corps qui transmet
la chaleur et celui qui la reçoit, comme c’est le cas de l’énergie rayonnée par le soleil. Ce sont
les photons, « particules » se déplaçant à la vitesse de la lumière, qui émettent de l’énergie
en fonction de leurs longueurs d’onde.
Modélisation à l’échelle d’un volume élémentaire
Le système à modéliser est décomposé en volumes élémentaires homogènes noté Vi .
Soit Ti la température du volume Vi , le bilan thermique de cette unité est décrit par l’application
de l’équation de conservation de l’énergie :
dTi
(ρC)i .Vi . = ∑φik (t) + Pi (t) (2.4)
dt k
avec :
(ρC)i : la capacité calorifique du matériau qui compose le volume Vi ,
Pi (t) : l’ensemble des sollicitations externes,
φik (t) : les flux échangés aux interfaces de Vi
Conditions aux limites
La résolution de l’équation de la chaleur à l’échelle du bâtiment nécessite de prendre en compte
les conditions aux limites du système analysé. Ces dernières se situent au niveau des murs, des
vitres et des planchers. Ils matérialisent l’interaction entre la surface et l’environnement sous
forme d’échanges convectifs et radiatifs. Parmi les sources de non linéarité, citons :
1. Pour les échanges convectifs de la surface extérieure de la paroi : le flux échangé est
φCE = hCE (TAE − TSE )
dépendant du coefficient d’échange convectif, de la température sèche de l’air extérieur et
de la température de la surface extérieure de la paroi. Avec hCE = a + bvn , fonction affine
de la puissance (non forcément entière) du vent.
2. Pour les échanges radiatifs : vient de la décomposition géométrique complexe de l’en-
soleillement (formules de Fresnel) au niveau des vitrages pour les échanges à courtes
longueurs d’ondes. En ce qui concerne les grandes longueurs d’ondes, le flux est égale-
ment non linéaire en fonction des températures (en T 4 , corps noirs).
Nous précisons que la discrétisation temporelle est souvent faite par la méthode des différences
finies (voir les travaux de Dautin [13] par exemple), mais peut également être obtenue par
intégration de l’équation 2.4.
19
2.1.2.3 Modèles réduits
La méthode de réduction modale Si le bâtiment est représenté sous forme d’état décrite
par l’équation 2.5, la réduction sous forme modale consiste à effectuer un changement de base
(base modale) pour diagonaliser la matrice d’évolution A (Dautin [13], Gao et al [19], etc.).
Xk+1 = AXk + BUk
(2.5)
Yk = CXk + DUk
Les valeurs propres les plus significatives sont gardées, cette opération est appelée troncature du
système dans sa base modale. L’intérêt de cette approche est la conservation des constantes de
temps du système initial dans le modèle réduit. Mais elle ne conserve pas les modes qui sont les
plus sensibles aux sollicitations et qui agissent le plus sur la sortie (notions de commandabilité
et d’observabilité).
Cette méthode, dite de Marshall ([42]), ne s’appuit que sur un critère temporel (constantes de
temps) pour effectuer la troncature. Cette dernière peut être améliorée par l’utilisation de la
méthode de Moore ([50]).
20
Le choix du nombre de résistances et capacités et leur configuration permet de restituer prio-
ritairement les basses fréquences du système et son régime permanent. L’analogie électrique
permet de comprendre la structure du modèle par une observation rapide du circuit équivalent.
Il est aisé de la modifier pour prendre en compte des spécificités du bâtiment (point d’injection
du flux solaire, couplage avec une autre zone thermique, etc.).
21
composante correspond au flux réfléchi par le sol : les flux direct et diffus sont en effet ré-
fléchis sur le sol de façon diffuse.
La proportion du flux diffus par rapport au direct est très variable. Pour des journées nua-
geuses, le direct peut être très faible et le diffus important.
ä La position relative du soleil par rapport à un point sur Terre, varie en fonction du temps.
C’est en partie cet aspect géométrique du rayonnement solaire qui introduit des termes dé-
pendant du temps dans les équations de la thermique du bâtiment.
Sortie P
Processus
-
Planification Entrée ε
d'expérience
+
Modèle
Sortie M
Méthode
Afin de cerner la réponse du processus sur une large bande fréquentielle, il est souvent question
d’utilisation de signaux riches fréquentiellement, dont le plus connu est la séquence binaire
pseudo aléatoire (SBPA), version numérique du signal blanc (contenant toutes les fréquences).
La seule limitation étant au niveau des conditions expérimentales.
Cette approche a fait l’objet des travaux de thèse de Richalet [58], pour différents types de
modèles : statiques et dynamiques d’ordre deux (réduits ou électriques équivalent). Les ensei-
gnements que l’on peut en tirer sont essentiellement :
ä L’introduction de paramètres physiques dans le vecteur des paramètres (constantes de temps,
GV,...) permet de simplifier l’initialisation de l’algorithme d’identification. En effet, le risque
que présente la plupart des algorithmes d’optimisation est celui de minimiseurs locaux, si
l’initialisation est trop éloignée de la réalité physique.
22
ä Les modèles statiques de signature énergétique nécessitent une période au moins hebdoma-
daire d’acquisition de données et une période d’expérimentation d’au moins dix semaines
alors que les modèles dynamiques permettent de travailler avec des périodes de l’ordre de la
semaine.
ä En l’absence de variabilité suffisante de l’ensoleillement, les contributions associées à la
sollicitation solaire sont systématiquement biaisées.
La détermination des caractéristiques dynamiques d’un procédé donné se fait dans le but de
concevoir un régulateur performant. Les conditions de fonctionnement normal sont en général
retenues et l’identification se fait selon le schéma suivant :
Processus réel
Données
Identification
Modèle
Commande
Régulateur
Les sollicitations ne font pas l’objet d’une planification particulière puisque le but est d’assurer
un fonctionnement optimal vis à vis des entrées réelles. Il en résulte que le modèle identifié est
adéquat pour ce type de scénario mais explique difficilement la dynamique du processus réel
pour un jeu de données différent. De même, le régulateur synthétisé est du même ordre que le
modèle (même niveau de complexité).
La procédure d’identification peut cependant être itérative (cas des commandes robustes), au-
quel cas il y a un couplage entre les étapes d’identification et de commande. En effet, une
erreur sur les paramètres du modèle induit une erreur sur les paramètres du régulateur, qui sont
par ailleurs des degrés de liberté du problème, puisque le critère optimisé porte sur le système
global régulé.
Parmi les travaux qui ont adopté cette approche en thermique du bâtiment, on peut citer :
ä Thèse de Alaoui (1992) [4] : la structure du modèle du bâtiment est de type fonction de
transfert d’ordre deux à trois entrées (température extérieure et flux solaire décomposé en
deux parties). L’identification est faite à deux échelles de temps : rapide et lente à partir de
filtrage passe haut et passe bas des données.
Chaque dynamique est identifiée dans son échelle de temps en conservant la continuité pour
23
les fréquences intermédiaires. La méthode a été appliquée à des données de synthèse et
les résultats sont plus satisfaisants que ceux obtenus aux moindres carrés classiques. Ces
résultats soulignent néanmoins la sensibilité de la méthode quant au choix de la fréquence
de coupure du filtre passe bas.
ä Thèse de Fraisse (1997) [18] : les travaux répondent à la question d’optimisation de la ges-
tion centralisée des scénarios d’intermittence. Le coût à minimiser est la consommation de
chauffage sans toute fois compromettre le degré de confort dans le bâtiment. La régulation
centralisée constitue une contrainte imposée par la réglementation thermique de 88 pour les
bâtiments à occupation discontinue. Cette catégorie regroupe typiquement les locaux ter-
tiaires dont le taux d’occupation ne dépasse pas 30% dans certains cas.
La commande floue est conçue, à partir d’un ensemble de quelques dizaines de règles
concernant les variables d’état, pour des cas test de bâtiments multi-zones. Ces bâtiments
sont modélisés sous forme de circuits électriques équivalents, d’ordre supérieur à 10 et dé-
passant la centaine de paramètres. La plus grande partie des paramètres est fixée à partir
d’une connaissance exhaustive du bâtiment et le calage du modèle par rapport au bâtiment
réel se fait au niveau de quelques variables comme l’inertie du bâtiment.
Bien que la solution développée soit plus intéressante que les relances à heure fixe ou en
fonction de la température extérieure, les performances du système de régulation par com-
mande floue montrent une grande sensibilité par rapport au nombre de règles floues utilisées.
Ceci reflète l’effet d’une incompatibilité entre la connaissance thermique servant à définir
les règles et le comportement thermique du bâtiment.
ä Thèse de Amstrong (2004)[3] : l’identification est faite pour des données enregistrées dans
des cellules test multi-zones, au pas de temps de 15 minutes. Plusieurs versions de modèles
de type ARX ont été envisagées selon le choix des entrées et sortie. En effet dans un premier
temps la puissance de chauffage/climatisation est une commande (ou entrée) alors que la
température intérieure est une observation (sortie), dans un deuxième temps c’est l’inverse
qui a été adopté et finalement, une combinaison de la température intérieure et de la consom-
mation de chauffage/climatisation est choisie comme observation.
Plusieurs ordres de modèle (de 1 à 5) ont été envisagés et les résultats soulignent surtout
la difficulté de choisir une température représentative pour un bâtiment multi-zone. Le troi-
sième type de modèle combinant la puissance et la température donne des résultats plus
satisfaisants que les deux premiersqui se sont revelés instables.
Le modèle identifié est ensuite utilisé pour la synthèse d’un régulateur de chauffage ou cli-
matisation avec un objectif double : maintien du confort et évitement des pics d’appel de
puissance. Cette démarche est menée dans un cadre de simulation et non sur des bâtiments
réels.
24
ä d’un régulateur du type Proportionnel Intégral (PI).
Néanmoins, la méthode d’identification proposée devrait être généralisable aux cas où le modèle
du bâtiment est un modèle d’ordre réduit autre que celui considéré dans la thèse.
Le choix du circuit équivalent R3C2, présentée sur la figure 2.6, pour modéliser le bâtiment est
motivé par le fait que, dans la catégorie des modèles de connaissance, il réalise un compromis
entre :
1. Les modèles de signature énergétique, généralement d’ordre 1 avec un gain et une constante
de temps interprétables physiquement. Ils ont été conçus pour décrire les variations gros-
sières de la consommation de puissance en fonction de l’écart entre la température exté-
rieure et une température intérieure fictive. Cette dernière peut par ailleurs correspondre à
un ensemble de bâtiment (à l’échelle du parc tertiaire bureaux par exemple) pour des pas
de temps élevés (journaliers, mensuels, annuels,...)
2. Les modèles électriques équivalents, plus détaillés, risquent d’être sur-paramétrés (beau-
coup de résistances pour décrire le détail des transferts dans les parois). Ces modèles,
développés sous forme de réseau, sont en général fonction de la géométrie du bâtiment.
Présentation et interprétation du modèle R3C2 Dans le modèle R3C2, présenté sur la fi-
gure 2.6, l’enveloppe du bâtiment est modélisée par l’ensemble (Ri , Ro , Cs ) (voir Marti et
Rignac[43], Jazouli et Barrois[31]). La résistance R f (en W /K) est une mise en parallèle de
trois résistances Rra , R p et Rle :
1 1 1 1
= + +
Rf Rra R p Rle
représentant respectivement les pertes dues au renouvellement d’air, les pertes hors renouvelle-
ment d’air (les deux types étant non déphasés) et la résistance des parois légères externes.
Les résistances Ri et Ro (en W /K) représentent les parois lourdes à forte inertie, le « pourcentage
Ri Ro » 5 donne plus ou moins d’importance à l’interaction avec l’intérieur ou l’extérieur avec le
déplacement du « potentiomètre » et le positionnement de la plus forte capacité du modèle, la
capacité de structure Cs .
Ri
5. Appellation utilisée par les thermiciens et qui désigne la quantité : Ri +Ro
25
Le flux solaire Qs est injecté directement au nœud de la température de structure Ts et la puis-
sance de chauffage Qres est injectée au nœud représentant la température intérieure Tint .
Le circuit électrique équivalent R3C2 est régi par les deux équations de premier ordre suivantes :
Synthèse
Les modèles développés pour décrire le comportement thermique du bâtiment se classent en
deux catégories selon que l’approche soit ascendante ou descendante. Dans le premier cas,
c’est le comportement statique qui à été d’abord la cible de la modélisation, menée avec des
pas de mesure mensuels, hebdomadaires ou journaliers. Ces modèles ont ensuite été enrichis
pour décrire le régime transitoire du bâtiment en incluant une constante de temps de quelques
dizaines d’heures, ce sont les modèles dits « de signature énergétique » du bâtiment étudiés avec
des pas de temps généralement horaires.
Dans le deuxième cas, le point de départ est une description détaillée du bâtiment à partir d’in-
formations très précises sur sa géométrie, la composition de ses parois et les données météoro-
logiques (température, flux solaire, vitesse du vent,...) au niveau de ses parois extérieures.
Des modèles d’ordre assez élevé (état à quelques centaines de composantes) sont obtenus à
partir de la résolution des équations de conservation de l’énergie sur des volumes de contrôle.
Selon la méthode de réduction adoptée, on peut soit obtenir des modèles d’ordres réduits (2
ou 3 généralement) de type boite noire (réduction de Moore, méthode modale par exemple)
soit des modèles dits « électriques équivalents » (si le modèle détaillé est mis sous forme de
circuit électrique analogue et que la réduction se fait en remplaçant des parties du circuit par
des conductances et capacités équivalentes).
26
Chapitre 3
27
Sollicitations Mesures
Planification Système réel
d'expérience
Connaissance
Ensemble de modèles
a priori
Critère
Invalidation
Oui
Non
La connaissance a priori : La connaissance a priori n’est pas une étape en soi, elle intervient
à plusieurs niveaux de la démarche d’identification à commencer par l’aspect physique du fonc-
tionnement du système réel. En effet, la connaissance des équations qui régissent la dynamique
du système permet de présélectionner une structure pour l’identification. Si de plus, le choix de
la séquence des données d’entrée est permis, la connaissance a priori des ordres de grandeur des
paramètres du modèle permet de mieux planifier l’expérimentation (fréquences d’excitation en
fonction des constantes de temps par exemple), afin d’optimiser la procédure d’identification.
Ensemble de modèles : Il existe deux classes de modèles, selon l’origine des équations dé-
crivant le comportement dynamique du système réel et l’abstraction des paramètres :
1. Les modèles physiques (ou de connaissance) : basés sur la décomposition du phénomène
en sous systèmes selon les lois de la physique. Les paramètres ont une interprétation
(caractéristiques des parois,. . . ). La mise en oeuvre de modèles de cette catégorie permet
d’exploiter des connaissances a priori fortes sur les valeurs des paramètres, mais elle
présente en général des non linéarités soit en l’état soit en les paramètres.
2. Les modèles comportementaux (ou de représentation) : le système est écrit sous forme de
relations purement mathématiques, qui expriment les sorties en fonction des entrées. Les
28
coefficients qui interviennent dans les équations n’ont pas de signification physique. Le
modèle ainsi obtenu est généralement appelé « boite noire ».
A(q-1)
u(k) y(k)
B(q-1)
S
F IGURE 3.2: Écriture en fonction de transfert
La relation entre l’entrée u et la sortie y du système S , dans ces conditions idéales d’absence
d’incertitude, est :
où les ai et bi sont les coefficients des polynômes A et B et q−1 l’opérateur retard (unitaire).
La partie déterministe de la structure est ainsi définie par le nombre de coefficients des deux
polynômes. Le vecteur θ des paramètres inconnus du modèle est égal à l’ensemble des coef-
ficients des polynômes dans le cas d’une modélisation en boite noire. Dans le cas général les
coefficients des polynômes A(q−1 ) et B(q−1 ) sont des fonctions des composantes deθ .
Il est parfois plus commode, lorsque la modélisation repose sur des équations de la physique,
de présenter le système S sous forme d’une représentation d’état :
X[k + 1] = A(θ )X[k] + B(θ )U[k]
Y [k] = C(θ )X[k] + D(θ )U[k]
29
Structure expression de la sortie inconnues
FIR y(t) = B(q−1 )u(t) + e(t) B
ARX A(q−1 )y(t) = B(q−1 )u(t) + e(t) A, B
ARMAX A(q−1 )y(t) = B(q−1 )u(t) +C(q−1 )e(t) A, B, C
B(q−1 )
OE y(t) = A(q−1 )
u(t) + e(t) A, B
−1
B(q ) C(q−1 )
BJ y(t) = −1
A(q )
u(t) + D(q−1 )
e(t) A, B, C, D
30
exemple), auquel cas l’estimateur θ̂ des moindres carrés est biaisé. Pour résoudre ce pro-
blème, on construit un « instrument » Z(t), corrélé avec les variables explicatives contenues
dans φ et non corrélé avec v. Le critère correspondant est dans ce cas :
N 2
J(θ ) = N1 ∑ Z(t) y(t) − φ T (t)θ
t=1
ä le maximum de vraisemblance : c’est une méthode statistique qui consiste à inférer les pa-
ramètres à partir de la série de mesures disponibles, y1:N . Si f (y1 , y2 , ..., yN |θ ) = f (y1:N |θ )
est la distribution des observations (ou mesures) par rapport aux paramètres, alors l’estima-
tion par maximum de vraisemblance consiste à trouver le jeu de paramètres qui permet de
N
maximiser cette fonction ou plus couramment son logarithme, J(θ ) = ∑ ln( f (yt |θ )) dans
t=1
le cas d’un bruit additif sur la sortie avec indépendances des yi .
Construction d’un modèle : C’est l’étape d’identification à proprement dite, elle consiste
à déterminer le modèle qui décrit le mieux le fonctionnement du système réel, parmi tous le
modèles dans la classe pré choisie. La construction du modèle est définie par le choix d’une
méthode ou algorithme d’identification. Elle fera l’objet des chapitres suivants.
Validation : C’est une étape importante qui permet de consolider les choix faits pour l’iden-
tification ou alors de remettre en cause certaines étapes dont :
ä le modèle (Walter et Pronzato [65]) : plusieurs aspects peuvent invalider le modèle comme
par exemple avoir des valeurs inadmissibles pour des paramètres physiques (résistance né-
gative,...), le caractère prédictif (une séquence d’entrée différente appliquée au modèle ne
permet pas d’approcher la sortie mesurée) ou la non robustesse (dépendance forte vis-à-vis
des conditions initiales ou vis-à-vis de perturbations faibles sur les entrées).
ä les hypothèses statistiques : certaines hypothèses sur les résidus analysés à travers des tracés
graphiques ou des tests statistiques de stationnarité, de gaussienneté ou de corrélation avec
les sollicitations, peuvent être remises en cause.
31
Dans un deuxième temps nous nous placerons dans le contexte restreint de notre étude c’est-
à-dire celui où les données collectées ne peuvent pas faire l’objet d’une planification. Nous
envisagerons donc les techniques qui permettent de quantifier la précision que l’on aura sur les
paramètres identifiés, compte tenu de la structure déjà choisie.
3.2 Identifiabilité
La question d’identifiabilité se pose au moment ou l’on choisit une structure. Elle s’effectue
généralement dans une étape préliminaire et permet de statuer, pour les observations obtenues
à partir du système réel et la structure envisagée, si le jeu de paramètres qui minimise le critère
est unique.
La préoccupation d’assurer au préalable que la structure M (θ ) envisagée ne donne pas lieu à
plusieurs solutions possibles, auquel cas la décision de choisir « la solution du problème » serait
délicate, a été étudiée sous plusieurs angles.
Dans le cadre du choix de la complexité de la structure, une première voie (Chappell [9],Walter
et Pronzato [65], Debus et al [14]) s’intéresse à la possibilité de différencier deux structures
M1 et M2 (figure 3.3) qui donneraient le même comportement de la sortie vis-à-vis des mêmes
sollicitations : c’est la notion de « discernabilité » (ou « distinguabilité ») par rapport à la sortie.
y1
M1 +-
u y
y2
M2 +-
Les structures présentées sont forme d’état sont ramenées à une écriture entrée/sortie avant de
faire l’objet d’une étude de discernabilité. Nous ne nous attarderons pas sur cette notion étant
donné que :
ä les problèmes d’identification sont généralement posés pour un modèle donné (comme c’est
le cas pour notre étude) ;
ä la « discernabilité » n’est ni une condition nécessaire ni suffisante pour l’unicité de la solu-
tion du problème d’estimation (Walter et al [64]).
Ainsi nous venons à la deuxième voie qui se focalise sur la propriété structurelle se rapportant à
l’unicité de la solution d’estimation : l’identifiabilité. Cette propriété a fait l’objet de plusieurs
études en fonction de la complexité de la structure envisagée. En effet, il n’existe pas de méthode
générale pour démontrer l’identifiabilité d’une structure non linéaire (la non linéarité peut être
relative à la dynamique et/ou la paramétrisation), mais un ensemble d’outils adaptés à chaque
situation et qui peuvent être classifiés en deux catégories :
32
1. algébriques : séries de Taylor (la sortie et ses dérivées successives), dérivées de Lie
(Walter et Pronzato [66]), méthode de changement de base (similarity transformation ap-
proach) (Chappell [9], Chappell et Gunn [10]), etc. Ces méthodes permettent de conclure
de manière certaine sur l’identifiabilité globale mais peuvent nécessiter des calculs com-
plexes.
2. numériques : permettent généralement de résoudre les problèmes complexes auxquels les
méthodes algébriques sont difficilement applicables. Les méthodes numériques apportent
généralement des réponses sur l’identifiabilité locale (i.e sur un voisinage d’une solution
et non sur tout le domaine des valeurs admissibles). Parmi ces méthodes nous citons les
approches basées sur l’analyse d’intervalles (Braems et al [7], Lagrange et al [36]) et le
calcul de noyau matriciel (Ollivier et Sedoglavic [54]).
Pour les besoins de notre étude, nous étudierons le formalisme de l’identifiabilité dans le cadre
particulier de structures de dynamiques linéaires, dans le domaine discret et présentées sous
forme de relations entrée/sortie.
3.2.1 Définition
Plusieurs définitions de l’identifiabilité existent selon l’objectif de l’approche d’identification.
En effet, soit le but est de garantir que l’estimateur θ̂N fourni par l’identification à partir d’une
séquence de mesures de longueur N converge vers le « vrai » θ0 qui est celui du système réel,
c’est la question de consistance de l’estimation qui se pose. La deuxième approche vise à ga-
rantir que l’estimateur obtenu est bien le minimum global du critère d’identification, c’est cette
approche d’unicité qui nous intéresse tout particulièrement. Les définitions présentées ici se
basent sur les travaux de Ljung et Gevers, le critère à minimiser est l’erreur de prédiction, mais
la démarche peut facilement être adaptée à d’autres critères.
Soit un modèle M (θ ) dont la relation entrée/sortie (Bazanella et al [6]) est décrite par :
où les fonctions de transferts peuvent être matricielles (cas de système multi entrées et/ou multi
sorties), q−1 est l’opérateur retard unitaire et e un bruit blanc centré (on remarque que toutes les
structures linéaires du tableau 3.1 peuvent se mettre sous cette forme). θ représente l’ensemble
des paramètres intervenant dans les coefficients des fonctions de transfert et qui sont l’objet de
la démarche d’identification. On notera que le modèle est complètement défini par la donnée du
−1 −1
couple G(q , θ ) H(q , θ ) .
La prédiction à un pas de la sortie en fonction des entrées et sorties est :
33
e
u G ++ y
Les transferts Wu et Wy peuvent être identifiés à partir de la structure dans le cas général présenté
ci-dessus : Wu (q−1 , θ ) = H −1 (q−1 , θ )G(q−1 , θ ) et Wy (q−1 , θ ) = I − H −1 (q−1 , θ )
Processus y
-
Entrée ε
Modèle y^ +
Identification
1 N T
θ̂N = arg min ∑ ε (t, θ )ε(t, θ )
θ N t=1
L’étude de l’identifiabilité a pour but de vérifier que le minimiseur du critère est unique sans
garantir que le vrai système admet la même structure que celle adoptée pour la démarche d’iden-
tification. Partant de cette constatation, l’identifiabilité locale d’une structure de modèle donnée
peut être quantifiée par le calcul du gramien d’identifiabilité :
34
ˆπ
∂W (e jw , θ ) ∂W H (e jw , θ )
Γ(θ ) = dw
∂θ ∂θ
−π
où Φz (w) est la densité spectrale de z. On remarque que la condition : ∀w, Φz (w) > 0, i.e. toutes
les fréquences sont présentes dans les signaux, est suffisante mais s’avère très restrictive.
Si l’expérience est informative et que la structure est globalement identifiable, ie chaque modèle
est représenté par un jeu unique de paramètres, la combinaison des deux propriétés permet
d’assurer une solution unique au problème d’identification. Il est néanmoins intéressant de noter
que :
ä Planifier une expérience se fait en général en vue d’identifier une structure donnée, or les en-
trées « designées » seront appliquées au système réel, a fortiori différent du modèle d’iden-
tification, ce qui implique une sortie différente de celle prévue lors de la planification de
l’expérimentation.
ä Il est toujours intéressant de comparer différentes structures, ce qui se fait généralement à
l’étape de validation. Pour cela il faudrait quantifier la richesse des signaux indépendamment
de la structure qui sera utilisée d’autant plus que dans certains cas, il n’y a pas ce degré de
liberté qu’est la planification d’expérience, auquel cas c’est plutôt le contenu des signaux
qui oriente le choix de la structure.
35
3.2.2.1 Contenu fréquentiel
Pour mener à bien l’estimation des paramètres d’un modèle, il est nécessaire de se demander
si les mesures acquises ou envisagées seront assez riches pour déterminer de manière assez
précise les paramètres. Cette préoccupation a notamment pour but de s’assurer que l’identité
du comportement du système physique d’une part et du modèle d’une autre part équivaut à
l’égalité des paramètres. Cette richesse se traduit par certaines propriétés du contenu fréquentiel
des signaux d’entrée en fonction de la structure et se quantifie par la notion de persistance.
T
Persistance : Un signal U(t − 1) , u(t − 1) u(t − 2) · · · u(t − n) est dit d’excitation
persistante à l’ordre n, si :
N
1
ä la limite de ru (τ) = lim ∑ u(t + τ)uT (t) existe et
N→∞ N t=1
ä sa matrice de covariance à l’ordre n :
ru (1) · · · ru (n − 1)
ru (0)
..
r (−1) ru (0) .
RUU = E U.U T = u .
.. ...
ru (1 − n) ··· ru (0)
est définie positive.
Pour illustrer la notion de persistance, on prend l’exemple d’un bruit blanc gaussien scalaire,
centré et de variance σ 2 . On a dans ce cas ru (τ) = σ 2 δτ0 , δ étant le symbole de kronecker (δ ji = 1
si i = j, 0 sinon). La matrice de covariance à l’ordre n est RUU = σ 2 In qui est toujours définie
positive. C’est cette propriété de persistance à tous les ordres qui explique l’utilisation fréquente
de signaux s’approchant de bruits blancs pour l’identification. On citera à titre d’exemple, la Sé-
quence Binaire Pseudo Aléatoire (voir figure 3.6), couramment utilisée en identification comme
« approximation » d’un bruit blanc discret, constituée d’impulsions rectangulaires (d’amplitude
a ou −a) modulées aléatoirement en longueur.
Richesse : un signal est suffisamment riche au degré n si sadensité spectrale est non nulle
pour au moins n fréquences distinctes de l’intervalle −π π .
La notion de richesse ne dépend que du signal et pas de la structure. Dans la réalité, il n’y a pas
de résultats sur le degré de richesse suffisant pour un problème d’identification donné mais on
peut déjà donner un degré minimal.En effet si le nombre de paramètres à déterminer est m, il
m+1
est nécessaire d’avoir au moins n = 2 fréquences, le signal le plus simple étant un mélange
de n sinusoïdes.
L’action du système à identifier sur chaque composante fréquentielle est caractérisée par le
déphasage et l’atténuation qu’elle induit. La comparaison de l’amplitude et de la phase, pour
chaque sinusoïde, à l’entrée et à la sortie du système permet de déterminer deux inconnues,
dans le cas où le système réel et le protocole expérimental sont idéaux (aucun bruit de mesure,
aucune perturbation parasite).
Dans la suite, on montre qu’à partir du calcul de la matrice d’information, on peut cerner les
limites de la précision maximale que l’on peut obtenir dans l’espace paramétrique. Ce critère
combine les caractéristiques d’identifiabilité de la structure avec le spectre des signaux d’en-
trée/sortie.
36
u(t)
-a
Si on suppose que le système réel S admet la même structure, c’est à dire qu’il existe θ0 tel
que S = M (θ0 ), alors l’erreur sur les paramètres converge vers une distribution gaussienne :
√
N(θ̂N − θ0 ) −→ N (0, I −1 (θ0 ))
N→∞
avec :
) ∂ ŷ(t|t−1,θ )
ä ψ(t, θ ) = − ∂ ε(t,θ
∂θ = ∂θ : jacobien de la prédiction du modèle
T
ä Σ = E e(t)e(t) est la matrice de covariance du bruit e.
Le calcul de la matrice d’information permet d’évaluer la borne inférieure de la variance de
l’erreur de convergence sur les paramètres. En effet, le cas présenté ici est optimal dans le
sens où le vrai système et le modèle d’identification admettent la même structure, mais aussi
parce qu’on se place dans le cas « idéal » où cette structure est globalement identifiable et que
les données utilisées sont suffisamment riches pour l’identification. Son calcul constitue un
indicateur sur la précision maximale que l’on peut avoir sur les paramètres identifiés compte
tenu du modèle et des données disponibles.
Pour mieux établir le lien avec les aspects d’identifiabilité et de richesse des signaux, nous nous
proposons d’exprimer la matrice d’information en fonction du spectre des données. En effet
ŷ(t|t − 1, θ ) = W (e jw , θ )z(t) d’où :
∂W (e jw , θ )
ψ(t, θ ) = z(t)
∂θ
37
z étant indépendant du vecteur des paramètre θ et par suite :
p ˆπ
1 ∂Wi (e jw , θ ) ∂W H (e jw , θ )
I(θ ) = ∑ Φz (w) i dw
i=1 2π ∂θ ∂θ
−π
y(t) = G1 (q−1 , θ )u1 (t) + G2 (q−1 , θ )u2 (t) + · · · + Gm (q−1 , θ )um (t) + H(q−1 , θ )e(t)
Plusieurs questions se posent quant à la démarche que l’on peut adopter pour estimer les coef-
ficients des m + 1 fonctions de transferts :
1. Vaut-il mieux procéder entrée par entrée et se ramener à m + 1 cas scalaires ou alors
exciter le système par toutes les sollicitations et le considérer dans sa globalité ?
2. Y a -t-il des structures particulières qui se prêtent plus au problème de l’identification
d’un système MISO ?
3. Que se passe t-il si des entrées sont corrélées entre elles ?
Ces questions ont fait l’objet de plusieurs études, et si la réponse à la première question a été
assez rapide la quantification de l’impact de la corrélation entre les données d’entrée s’avère
beaucoup plus subtile. En effet les hypothèses sous-jacentes aux études menées écartent le cas
de corrélation entre les sollicitations et sont :
ä les entrées u1 , · · · , um sont statistiquement indépendantes entre elles,
ä chaque signal ui est d’excitation persistante d’ordre ni , si ni est la somme des degrés des
polynômes constituant la fonction de transfert Gi .
Sous ces hypothèses, les principaux résultats établis sont :
1. l’ajout d’une entrée supplémentaire um+1 , vérifiant les hypothèses précédentes, relative
à une fonction de transfert Gm+1 améliore la précision de tous les paramètres communs
entre Gm+1 et les autres fonctions de transferts (voir Gevers et al[24]). Cette démons-
tration consolide le choix de l’excitation simultanée de toutes les entrées, plutôt que de
considérer plusieurs identifications de sous systèmes mono entrée mono sortie.
2. l’expression asymptotique de la variance de l’estimée du bruit :
n Φv (w)
var(ĤN (e jw )) w
N |H(e jw )|2
38
r u y
+- Régulateur Système
Capteur
où v est le bruit qui s’additionne à la sortie du système (voir figure 3.8) et n l’ordre de H
a été établie depuis longtemps par Ljung [40]. L’expression exacte de cette variance a été
calculée dans la cas des structures BJ et OE par Ninness et Hjalmarsson [53] qui montrent
son indépendance des entrées ui . Ces résultats ne sont cependant pas valables pour les
structures ARX et ARMAX pour lesquelles Gevers et al [24] montrent que chaque en-
trée contribue à diminuer cette variance ainsi qu’à augmenter la précision des transferts
associés aux autres entrées.
avec cette fois : u(t) = r(t) − F(q−1 )y(t), en notant S(q−1 , θ ) = 1+G(q−11)F(q−1 ) la sortie y s’ex-
prime en fonction des entrées indépendantes r et v de la manière suivante :
v étant le bruit additif sur la sortie du système en boucle ouverte : v = H(q−1 )e(t) comme
l’illustre la figure 3.8. La dépendance des fonctions des transferts en les paramètres a été omise
pour alléger les notations.
De même la commande u s’exprime en fonction de r et v : u(t) = S(q−1 )r(t)−F(q−1 )S(q−1 )v(t),
on en déduit son spectre :
2 2 2
Φu (w) = S(eiw ) Φr (w) + F(eiw ) S(eiw ) Φv (w)
Il est clair que dans le cas où l’action du régulateur n’est pas nulle, il y a corrélation entre
l’entrée et le bruit qui n’est pas mesurable (ces résultats ont été établis par Akaike [1]). Cette
39
e
v
r u + y
+ G +
-
Système en BO
corrélation résulte typiquement en un biais sur l’estimation, celui-ci a fait l’objet de plusieurs
études qui visent à le transformer pour servir le but de synthèse de régulateur (voir Van Den
Hof et Schrama [30] par exemple).
L’identification des paramètres inconnus en boucle fermée passe par le choix du système ou sous
système à identifier. Plusieurs cas de figures existent en fonction de l’information disponible sur
le régulateur F :
1. Aucune connaissance n’est disponible pour le régulateur, dans ce cas la rétroaction est
ignorée et une approche boucle ouverte ne prenant pas en compte le signal de référence
r, même s’il est connu, peut être menée. On peut par exemple appliquer la méthode de
l’erreur de prédiction au système asservi en ne considérant que u et y, et plus généralement
toutes les méthodes d’identification en boucle ouverte à condition de bien choisir les
séquences de données ou de les filtrer le cas échéant (voir Landau [37]). C’est l’approche
directe.
2. Le régulateur ainsi que le signal de commande sont connus : on identifie d’abord la boucle
fermée en prenant le système global ayant pour entrée r et pour sortie y, plusieurs algo-
rithmes spécifiques de l’identification en boucle fermée peuvent être utilisés (voir à titre
d’exemple Landau [37]). Dans un deuxième temps, on détermine les paramètres du sys-
tème en boucle ouverte en se basant sur les estimations de la boucle fermée. Il s’agit dans
ce cas de l’approche indirecte.
3. Le régulateur est inconnu mais de structure connue, dans ce cas on prend u et y comme
sorties d’un système d’entrée r. Plusieurs méthodes permettent ensuite de retrouver les
paramètres de la boucle ouverte. C’est l’approche jointe entrée-sortie.
L’étude détaillée de ces approches dans le contexte de minimisation de l’erreur de prédiction a
fait l’objet des travaux de thèse de Forssell [17]. Les principaux enseignements qui en découlent
sont :
ä Si le but de la régulation est de stabiliser le système, c’est à dire de diminuer sa sensibilité par
rapport aux fréquences présentes dans les perturbations, cette propriété contribue fortement
à la dégradation de l’informativité des données surtout en ces fréquences.
ä Plusieurs méthodes utilisées pour l’identification de systèmes opérant en boucle ouverte
comme les méthodes spectrales, la variable instrumentale et les méthodes à erreur de sortie
donnent des résultats médiocres si elles sont appliquées à l’approche directe avec un modèle
de bruit (par rapport au système réel) incorrect.
40
ä L’approche indirecte est performante mais risque de donner de très mauvais résultats (plus
mauvais que ceux de l’approche directe) si le modèle du régulateur, qui est supposé connu,
n’est pas bon.
Synthèse
Dans ce chapitre, nous avons présenté la procédure d’identification en mettant le trait sur les
étapes préliminaires de choix de la séquence des données et de la structure du modèle à identi-
fier.
Les paramètres de ce modèle, de structure fixée, sont estimés par rapport à la capacité de celui-ci
à prédire le comportement du système réel. Mais pour que la solution du problème d’estima-
tion correspondant soit unique il faut vérifier deux propriétés : l’identifiabilité de la structure et
l’informativité des données. Pour statuer sur le caractère identifiable/non identifiable de la struc-
ture considérée, nous avons introduit les outils théoriques d’étude de l’identifiabilité. Quant à
l’adéquation entre la complexité de la structure et de la richesse des données, qui conditionne
la qualité de l’estimation, nous avons présentés quelques outils qui permettent de quantifier la
précision sur les paramètres estimés.
Par ailleurs, la complexité du problème de l’identification est augmentée quand le système est
multi-entrées ou en boucle fermée. C’est à ce titre que nous avons présenté une synthèse de
travaux antérieurs décrivant les considérations à prendre en compte pour identifier ces familles
de structures.
41
Chapitre 4
Structure Critère
Détermination
du modèle
Invalidation
43
1. Une structure incluant la dynamique du système et celle du bruit qui traduit les incerti-
tudes
2. Un critère de sélection qui permet de choisir, dans la classe déjà fixée par le choix de
la structure, le modèle dont la dynamique explique le mieux les informations contenues
dans les signaux.
3. Un algorithme qui permet de minimiser ce critère pour obtenir un jeu de valeurs numé-
riques optimales des paramètres.
Ainsi pour mener une estimation paramétrique dans le cas de la thermique du bâtiment, nous
nous proposons d’abord de choisir une représentation de la structure (entrée, sortie, état) adaptée
au système global « bâtiment + régulateur ». Cette structure est ensuite complétée par le choix
de la dynamique du bruit. Enfin, nous décrirons brièvement l’algorithme utilisé pour minimiser
le critère sachant que nous nous plaçons dans le cadre de minimisation de l’erreur de prédiction.
uobs yobs
Capteur 1 Capteur 2
r r - yobs u y
+- Régulateur Bâtiment BO
44
{}
Text
AI Bâtiment BO
+ Qch
Qs
E Régulateur
F IGURE 4.3: Système global « batiment + régulateur » ayant comme entrées : la température
extérieure Text , les apports internes AI, le flux solaire Qs et la consigne de température E et
comme sortie la courbe de charge chauffage Qch
E +- Régulateur Qch
++ AI
F IGURE 4.4: Modèle du bâtiment en boucle fermée, les températures de consigne E et extérieure
Text , le flux solaire Qs et les apports internes AI sont des entrées du modèle global, Qch en est
la sortie.
de transfert reliant les sollicitations (température extérieure, apports internes, flux solaire et
température de consigne) à la sortie courbe de charge chauffage comme l’illustre la figure 4.3
Cette limitation du point de vue mesure (indisponibilité de la température intérieure) induit une
grande complexité de la paramétrisation du système global. Même si ce dernier est linéaire par
rapport à l’état et d’ordre faible, toute la question se pose vis-à-vis de son identifiabilité et de
la sensibilité de la courbe de charge chauffage vis-à-vis de la variation de chaque paramètre du
modèle.
Nous rappelons que le modèle choisi pour le bâtiment est le circuit électrique équivalent R3C2
(illustré sur la figure 4.5), celui du régulateur est de type PI, interagissant selon le schéma de la
figure 4.4.
Modèle entrée/sortie du système global La formulation de la boucle fermée est basée sur les
équations du modèle R3C2 présenté précédemment et la prise en compte de la rétroaction, qui
décrit le fonctionnement du régulateur PI :
h i
dTint 1
dt = Cres − R1f + R1i Tint + R1i Ts + R1f Text + Qres
h i
dTs 1 1 1 1 1 (4.1)
dt = Cs Ri Tint − Ro + Ri Ts + Ro Text + Qs
´
Qch = K E − Tint + τ1 (E − Tint )
La troisième équation du système 4.1 décrit la dynamique du PI, ce dernier est caractérisé par
son gain K et son temps d’intégration τ. En choisissant comme composantes du vecteur X
45
F IGURE 4.5: Schéma électrique équivalent R3C2
´
d’état : Tint , Ts et d = τ1 (E − Tint )(état de l’intégrateur), le système précédent se met sous
forme d’état compatible avec la modélisation de la figure 4.3 :
Ẋ(t) = AX(t) + BU(t)
Qch (t) = CX(t) + DU(t)
Ce système est ensuite discrétisé par la méthode d’Euler explicite, ce qui nous place dans le
cadre de systèmes discrets pour la suite du chapitre (les matrices de la représentation d’état
discrète seront notées Ad , Bd , C, D ; ces deux dernières sont invariantes par passage continu/-
discret).
46
Cette représentation, bien que fonction de deux bruits différents, est équivalente à :
Xk+1 = Ad Xk + Bd Uk + Mek
(4.2)
yk = CXk + DUk + ek
où M 1 est le gain de Kalman associé au système stationnaire (voir Stoica et Soderstom pour la
démonstration qui utilise le théorème de la factorisation spectrale). M est fonction des matrices
de covariance des bruits v et e. La représentation 4.2 est appelée forme innovation.
Par rapport aux notations du chapitre précédent, les fonctions de transfert associées à la structure
entrée/sortie s’identifient de la manière suivante :
Le choix d’isoler le terme H0∗ (q−1 )ek est dû au fait qu’il est fonction de l’erreur à l’instant k − 1
et aux instants antérieurs. Il est donc parfaitement connu à l’instant k. Étant donné que ek est
décorrelé par rapport aux données et indépendant des échantillons passés, la variance de l’erreur
de prédiction est minimale lorsque :
47
ŷk = G0 (q−1 )uk + H0∗ (q−1 )(yk − ŷk )
= H0−1 (q−1 )G0 (q−1 )uk + (I − H0−1 (q−1 ))yk
En remplaçant H0 et G0 par leurs expressions de la représentation d’état, on montre que le
modèle récursif d’estimation est :
X̂k+1 = (Ad − MC) X̂k + Bd Uk + Myk
(4.3)
ŷk = CX̂k + DUk
Comme on s’est placé dans un cadre de minimisation de l’erreur de prédiction, le critère à
minimiser dans ce cas est :
N
1 2
J(θ ) = N ∑ kyk − ŷk k
k=1
N 2
= 1
N ∑ H0−1 (q−1 )[yk − G0 (q−1 )uk ]
k=1
où y est Qch et que u est composé de : Text, Qs, AI et E.
48
Algorithme 4.1 Algorithme de minimisation de l’erreur de prédiction pour le système sous
forme innovation
1. Initialisation, θ̂0
2. Changement de base équilibrée : θ̂0 → θ̂0e
h 2
i
e
3. Résolution de θ̂k = arg min J(θ ) + δ2 θ − θ̂k−1 par Gauss-Newton (voir annexe B)
θ
4. θ̂k → θ̂ke : réalisation équilibrée
5. Reprise à l’étape 3
Cres 5
Cs 8
et K = 2
ä Paramètres du modèle : Rf = 5
τi 4
Ri 2.5
Ro 0.5
16
re f
ä Période d’échantillonnage et état initial : ∆t = 1 s, X0 = 15
5
ä Variances des bruits d’état et de mesure : Q = diag( 0.16 0.09 0.09 ), R = 1
Comme les résistances, capacités et temps d’intégration associés au modèle d’état global du
bâtiment interviennent par inverse, nous avons effectué le changement de variable décrit dans
le tableau 4.1 :
zres = 1/Cres
zs = 1/Cs
zτ = 1/τi
z f = 1/R f
zi = 1/Ri
zo = 1/Ro
zK = K
Tableau 4.1: Changement de variable du modèle global
49
F IGURE 4.6: Écart entre courbe de charge de référence et simulée, pour une perturbation de
1 1
référence ref_bruit : ΣX2 0 = diag([0.2; 0.1; 0.2]), Q 2 = diag([0.4; 0.3; 0.3]), R = 1
L’identifiabilité locale se rapporte à l’étude des variations de la sortie lorsque les paramètres va-
rient dans un voisinage infinitésimal autour des paramètres de référence (ou point d’équilibre).
Les directions paramétriques les plus influentes sont celles dont les composantes correspon-
dantes du gradient du gramien d’identifiabilité sont les plus élevées.
L’étude de l’identifiabilité locale est menée dans le cadre du modèle à erreur de sortie (H = I
ou M = 0). Dans cas, le gramien d’identifiabilité vaut :
ˆπ
∂ G0 (e jw∆t , θ ) ∂ GH
0 (e
jw∆t , θ )
Γ(θ ) = dw
∂θ ∂θ
−π
w∆t
La période d’échantillonnage est ∆t , soit en fréquence réduite ν = 2π le gramien est défini
pour − 21 ≤ ν ≤ 12 .
Le gramien d’identifiabilité a été envisagé pour un jeu de paramètres θ type, qui peut servir
comme initialisation pour l’algorithme d’identification. Pour calculer Γ, on prend en considé-
ration le fait que la fonction à intégrer est paire en ν donc on restreint l’intervalle d’intégration
50
F IGURE 4.7: Écart quadratique entre courbe de charge de référence et simulée, pour une per-
turbation dans le sens du gradient, grad : gradient normalisé
à 0 , 12 qui a été discrétisé en 100 points. L’integrale est approchée par la méthode des tra-
pèzes 4 .
Le gramien d’identifiabilité Γ est une matrice hermitienne. Pour montrer l’identifiabilité locale
du problème, ilsuffit de vérifer que Γ est inversible, pour cela, ses valeurs propres sont présen-
tées ci-dessous :
-2.6352e-010
-1.6433e-010
7.5570e-011
3.0420e-010
3.8450e-001
1.7643e+000
1.4374e+001
Il est clair que ce gramien n’est pas inversible numériquement, seules trois directions (sur les
sept) peuvent être déterminées avec précision. Pour mettre en évidence les paramètres les plus
impactant sur ces directions, nous présentons les vecteurs propres associés.
51
600
dzres
dzs
500
dzτ
dzf
400
dzi
dzo
300
dzk
δY/δθ
200
100
−100
0 10 20 30 40 50 60 70 80
Temps (h)
F IGURE 4.8: Sensibilité Qch par rapport à une variation de chacun des paramètres
Ainsi, il en ressort que le paramètre dont la précision sera la meilleure est celui dont la contribu-
tion au dernier vecteur propre est la plus considérable : soit l’inverse de la capacité de structure
Cs .
Afin de quantifier la sensibilité de Qch par rapport à une variation infinitésimale de chacun des
)
paramètres, nous nous proposons d’évaluer ∂ ŷ(t,θ
∂θ où y est Qch . Pour cela on pose X = ∂∂ X̂θ et
Y = ∂∂θŷ , la différentiation de 4.3, en prenant un terme d’innovation nul (M = 0), par rapport
aux paramètres donne :
X̂k+1 = Ad X̂k + Bd Uk
X = Ad Xk + ∂∂Aθd X̂k + ∂∂Bθd Uk (4.4)
k+1
Yk = CXk + ∂ θ X̂k + ∂ θ Uk
∂C ∂D
52
D’une autre part, les variations de l’inductance zi et du temps d’intégration du régulateur zτ
semblent impacter peu Qch . Ceci peut être la conséquence du fait que le terme proportionnel est
prépondérant dans l’action du régulateur d’un coté et d’un autre coté, une fois chargées (valeur
conséquente pour la première période), les parois légères du bâtiment ne conduisent plus la
chaleur.
20
19
0 10 20 30 40 50 60 70
8
Text (°C)
6
4
0 10 20 30 40 50 60 70
Qs (kW)
0.5
0
0 10 20 30 40 50 60 70
0.6
0.5
Qch (kW)
0.4
0.3
0.2
0.1
0 10 20 30 40 50 60 70
Temps (h)
Le bâtiment modélisé est constitué de cinq zones thermiques soumises à la même consigne de
température (intermittence 19/20°C) et au même type de système de chauffage, l’une des zones
présente une saturation de la puissance de chauffage (courbe en magenta sur la figure 4.9).
Le système équivalent considéré est supposé avoir une température intérieure, non disponible,
obtenue comme moyenne des cinq zones et une consommation égale à la somme des consom-
mations de toutes les zones. Une analyse rapide des données fait ressortir une périodicité jour-
nalière des signaux de sollicitations comme l’illustre la figure 4.10.
Pour mieux faire ressortir cette périodicité, on a calculé les densités spectrales de puissance
associées aux signaux, illustrées dans la figure 4.11. En effet, la période de 1 jour semble pré-
pondérante dans les spectres, qui ont par ailleurs un contenu fréquentiel très pauvre en hautes
fréquences. Ces fréquences sont celles dont le module du spectre est inférieur à 0 dB, soit pour
les périodes temporelles inférieures à 4 heures, à l’exception du flux solaire.
53
8
Text(°C)
6
4
0 0.5 1 1.5 2 2.5 3
1
Qs (kW)
0.5 Variations journalières
0
0 0.5 1 1.5 2 2.5 3
20
E (°C)
19
F IGURE 4.10: Données du bâtiment considéré comme une zone thermique unique
60 60
Φ
40 Text
20 ΦQs
40 0
−20 ΦE
Densité Spectrale de Puissance (dB)
−40
20 0 2 4 6 8 Φ
Qch
−20
−40
−60
−80
54
0.25
Qch mesurée
Qch(0.1KW)
0.15
0.1
0.05
10 20 30 40 50 60 70
Temps (h)
θ0 θf : M = 0 θ f : M 6= 0
7 . 4 0 4 6 e −007 9 . 2 1 7 8 e −007 1 . 0 3 6 5 e −006
6 . 4 9 3 7 e −012 7 . 1 6 7 5 e −008 7 . 1 1 5 3 e −008
0.0032256 0.0030396 0.0059622
1172.6 19.979 29.193
2496.9 400.38 446.73
857.27 81.191 67.345
856.75 87.787 52.874
Tableau 4.2: Résultat d’identification des paramètres avec modèle de bruit estimé
55
entrées/sortie adopté pour l’identification est :
Text [k]
Qch [k] = G θ , q−1 Qs [k] + e [k]
E [k]
où e : erreur sur la sortie est également le résidu d’identification puisque la courbe de charge
estimée par le modèle est :
Text [k]
Q̂ch [k] = G θ̂ , q−1 Qs [k]
E [k]
Le résidu e est considéré par hypothèse comme une réalisation d’un processus aléatoire gaussien
de moyenne nulle et de variance donnée Σ.
Ces tests ont pour but de vérifier si les hypothèses faites sur les résidus sont respectées, et parmi
ces hypothèses, on s’intéresse plus particulièrement à :
1. Caractère gaussien
0.999
0.997
50
0.99
0.98
0.95
40 0.90
0.75
Probabilité
30
0.50
0.25
20
0.10
0.05
0.02
10 0.01
0.003
0.001
0
−0.01 −0.005 0 0.005 0.01 0.015 0.02 −0.01 −0.005 0 0.005 0.01 0.015 0.02
erreur de prédiction: e erreur de prédiction: e
Synthèse
L’algorithme d’identification a été appliqué à un jeu de données provenant d’un logiciel de
simulation du bâtiment où l’on modélise, sous forme détaillée, un bâtiment avec plusieurs zones
thermiques, dont une avec saturation de la puissance de chauffage.
56
Autocorrélation de ε
1.2
0.8
0.6
φεε(τ)
0.4
0.2
−0.2
0 5 10 15 20 25
τ
0.3
0.2
0.1
φεT (τ)
ext
−0.1
−0.2
−0.3
−0.4
−25 −20 −15 −10 −5 0 5 10 15 20 25
τ
L’évaluation de l’approche d’identification est appuyée par une étude structurelle du modèle
global à travers la quantification de son identifiabilité et la caractérisation de l’informativité des
signaux de sollicitation. Les résultats obtenus présentent une grande sensibilité par rapport à
l’initialisation de l’algorithme, ce qui est révélateur du caractère mal posé du problème.
57
Intercorrélation entre ε et Qs
0.4
0.2
φεQ (τ)
0
s
−0.2
−0.4
−25 −20 −15 −10 −5 0 5 10 15 20 25
τ
Intercorrélation entre ε et E
0.4
0.2
φεE(τ)
−0.2
−0.4
−25 −20 −15 −10 −5 0 5 10 15 20 25
τ
58
Chapitre 5
Introduction
Les différents tests de sensibilité menés au chapitre précédent ont montré que le caractère mal
posé du problème risque d’amener à des valeurs des paramètres identifiés non pertinentes. De-
vant ce constat la méthode d’identification doit être enrichie pour ne pas fournir seulement une
valeur des paramètres, et peut être une covariance, mais toute l’information sur leurs distribu-
tions apportée par les mesures et les connaissance ou choix a priori. Pour résoudre ce problème
nous utiliserons la simulation bayésienne.
L’objectif que nous proposons dans cette partie du mémoire est donc d’obtenir les distributions
a posteriori des paramètres du modèle R3C2 asservi en température par un régulateur Propor-
tionnel Intégral ainsi que des variables d’état cachés (températures intérieure et de structure).
La présentation du cadre bayésien et des méthodes classiques de simulation stochastique des lois
a posteriori feront l’objet de la première partie de ce chapitre. La deuxième partie sera consacrée
à l’application de l’approche bayésienne pour le bâtiment et la formulation de l’échantillonneur
de Gibbs qui s’en suit.
59
Dans cette approche d’inversion, répondre à la question d’identification du bâtiment consiste
à inférer l’entrée inconnue θ (paramètres du modèle dans notre cas) d’un système de paramé-
trisation H connue en fonction de θ et dont y est une observation, comme l’illustre la figure
5.1. L’inférence, qui contrairement à la modélisation, permet d’aller du domaine du tangible
(données, mesures) à un domaine plus abstrait présente deux difficultés au niveau de :
1. l’exploitation de l’information disponible dans les données et
2. la mise en forme du problème en vue de sa résolution mathématique.
y H(Θ, U) Θ
La structure H inclut les modèles du bâtiment et du régulateur de température, soit sept para-
mètres inconnus. Par ailleurs, la détermination des composantes cachées de l’état (températures
intérieure et de structure) aide à analyser le comportement thermique du bâtiment. En effet, la
comparaison entre la température de consigne -supposée connue- et la température intérieure
ainsi que la reconstitution de la température de structure constituent un moyen de validation
du modèle inféré d’un point de vue thermique. L’importance accordée à ces deux composantes
découle du sens physique direct qu’elles ont.
Compte tenu de ces considérations, la structure à inférer sera mise sous forme de représentation
d’état stochastique :
XT [k + 1] = A [θ ] XT [k] + B [θ ]U [k] +V [k]
Y [k] = C [θ ] XT [k] + D [θ ]U [k] +W [k]
La mesure (ou observation) Y est la courbe de charge de chauffage Qch et les sollicitations
connues regroupées dans le vecteur U sont respectivement : la température extérieure, le flux
solaire, les apports internes (calorifiques autres que le chauffage) et la température de consigne
(souhaitée par l’occupant). Les bruits d’état V et de mesure W traduisent respectivement l’in-
certitude par rapport au modèle choisi (ou dynamique) et celle par rapport à la précision des
mesures (capteur, acquisition des données,...).
La modélisation du problème pour l’inférence est illustrée sur la figure 5.2 :
où le vecteur des inconnues X regroupe les paramètres Xθ (c’est à dire θ ) et les séries de tempé-
ratures XT . Les informations a priori disponibles pour ces inconnues sont :
60
ä une table de valeurs types des résistances et capacités du modèle R3C2 pour des bâtiments
différents. Cette table est élaborée à partir d’études thermiques de bâtiments menées à EDF
(voir la fin de l’annexe A). Cette table permet de donner des ordres de grandeurs pour les
différents paramètres mais n’est pas suffisamment fiable pour suggérer l’utilisation de lois a
priori discrètes.
ä pour les bâtiments tertiaires (bureaux), envisagés ici, il est généralement possible de disposer
de document décrivant le système de chauffage. Par exemple, pour le bâtiment LGEP1 du
site des Renardières, il s’agit d’un ensemble de régulateurs PI de paramètres connus.
ä en ce qui concerne les températures, les profils possibles se situent généralement dans un
intervalle connu. Par exemple, dans le cas d’intermittence, la température intérieure varie
généralement entre 18 et 23 degrés.
L’exploitation de ces informations a priori sur les inconnues est une étape clé de l’inversion
bayésienne, elle constitue le point de départ pour l’apprentissage statistique des données. La
section suivante sera dédiée à la présentation dans un cadre général, celui du schéma 5.2, de
l’approche bayésienne.
61
2. la vraisemblance des données mesurées π(Y |X) : elle quantifie la probabilité théorique
des observations à partir des entrées étant supposé connues les inconnues X. Cette loi de
probabilité est en général obtenue en faisant l’hypothèse de bruit additif indépendant de
X:
Y = H (X) + E
E étant le modèle de bruit -ici linéaire en V et W - dont la probabilité de distribution est
supposée connue πbruit . On en déduit la vraisemblance :
La loi conjointe de X et Y permet de mixer ces deux sources d’information grâce à la relation
de Bayes :
π(Y |X)π(X)
π(X|Y ) = ´
π(X)π(Y |X)dX
Cette probabilité qui combine connaissance a priori des inconnues et les informations contenues
dans la mesure est la distribution a posteriori de X. Elle peut servir à estimer plusieurs grandeurs
caractéristiques de X a posteriori : moments, quantiles,. . . En effet, toute fonction g de X peut
être estimée a posteriori de la manière suivante :
´
g(X)π(Y |X)π(X)dX
E [g(X)] = ´
π(X)π(Y |X)dX
Plusieurs difficultés peuvent apparaître à ce niveau :
1. le calcul explicite de la loi a posteriori peut être impossible,
2. la fonction g est compliquée
3. même si ces lois sont connues, lorsque l’intégrale n’a pas une forme analytique, le calcul
numérique de E [g(X)] peut nécessiter une puissance de calcul considérable surtout en
grande dimension.
Ces difficultés d’ordre pratique sont à l’origine du développement de techniques de simulation
stochastique qui constituent un moyen d’approximation basées sur la méthode de Monte Carlo.
Ces approches feront l’objet de la section suivante où on abordera particulièrement les méthodes
à noyaux de Metropolis Hastings qui seront utilisées par la suite.
62
5.1.2 La simulation stochastique
On entend par simulation dans le contexte de notre étude la mise en évidence de la variabilité
résiduelle engendrée par π (Y ). Cela se fait par la génération de scénarios à partir des lois
conjointes définies au 5.1.
Nous nous intéressons ici plus particulièrement à la simulation de la loi a posteriori. La généra-
tion d’échantillons directement à partir de la loi a posteriori n’est pas facile sauf dans certains
cas standards. Par exemple, dans le cas où la loi est scalaire, on peut être amené à utiliser des
techniques directes 1 comme le changement de variable à partir d’une loi simple dont l’échan-
tillonnage est maîtrisé (loi uniforme typiquement).
Dans les cas multidimensionnels et/ou non standards, on fait appel à l’échantillonnage par
chaîne de Markov. L’idée principale est de générer une suite de points distribués asymptoti-
quement selon la loi a posteriori. Pour générer ces échantillons, on fait appel aux techniques
MCMC, faciles, efficaces et basées sur l’utilisation de la méthode de Monte Carlo.
Ayant à disposition cette suite {xi , i = 1 · · · n} qui admet pour limite ergodique la loi cible (loi
a posteriori dans le cas envisagé ici), la méthode de Monte Carlo consiste à approcher la valeur
exacte de l’intégrale en évaluant la fonction à intégrer en ces points :
ˆ
1 n
g(X)π(X|Y )dX = E [g(X)] ≈ ∑ g(xi )
n i=1
Lorsque les points xi forment une chaîne ergodique, la loi des grands nombres assure la conver-
gence de la somme vers la valeur de l’intégrale quand le nombre d’échantillons augmente.
Pour comprendre les techniques de MCMC qui fournissent la suite xi , on présente la définition
de chaîne de Markov ainsi que les propriétés nécessaires sans rentrer dans le détail de la théorie
(on pourra consulter Feller [16] à ce titre).
Il s’agit d’un processus stochastique {x0 , · · · , xn } à temps discret avec la propriété que l’infor-
mation induite par la connaissance de toutes les variables antérieures à t est résumée à l’instant
t − 1 par xt−1 (t − k à t − 1 pour une chaîne d’ordre k, on considère ici les chaînes d’ordre 1), ce
qui se traduit en terme de densité de probabilité de la manière suivante :
P(xt |x0 , x1 , · · · , xt−1 ) = P(xt |xt−1 )
On définit le noyau de transition de la chaîne, Kt (.|.) = P(xt |xt−1 ). Lorsque ce noyau est indé-
pendant de l’indice t la chaîne est dite homogène.
Pour une chaîne homogène, le noyau de transition peut être construit de manière à laisser in-
variantes certaines lois. On entend par invariance d’une densité π ∗ par K (.|.) le fait que :
π ∗ = K π ∗ , ce qui implique que si xt−1 ∼ π ∗ , il en sera de même pour son successeur xt .
Partant de cette propriété, résoudre le problème de simulation stochastique par chaînes de Mar-
kov revient à construire un noyau de transition qui laisse la loi a posteriori invariante, tout en
garantissant la convergence de la chaîne ainsi construite vers la loi a posteriori.
1. D’autres techniques existent comme les méthodes de mélanges, la méthode d’acceptation-rejet ou les mé-
thodes générales pour les densités log-concaves.
63
D’autres propriétés plus techniques comme :
ä l’apériodicité : une chaine est apériodique si elle n’est pas cyclique, Robert [60]) et,
ä l’irréductibilité par rapport à π ∗ : indépendamment du point initial, une chaîne irréductible
permet d’explorer tout l’espace des points (ayant une probabilité non nulle pour π ∗ )
sont nécessaires pour assurer la convergence en loi de la somme vers l’intégrale dans le cas
fini 2 .
Les techniques MCMC sont apparues en 1953 (Metropolis et al [49]) pour la physique statis-
tique et en 1990 (Gelfand et Smith [20]) dans la littérature statistique grâce au développement
des moyens informatiques (on renvoie le lecteur à Gilks et al [26] par exemple pour la mise
en oeuvre pratique des méthodes MCMC). Elles utilisent des noyaux de transition réversibles
par rapport à la loi cible (c’est à dire : π ∗ (x)K (y|x) = π ∗ (y)K (x|y)), propriété plus commode
impliquant la stationnarité. Parmi ces méthodes, les algorithmes de Metropolis-Hastings et de
Gibbs sont les plus connus. Ces algorithmes paraissent adaptés au problème de caractérisation
du bâtiment respectivement dans le cas de modélisation en temps continu et en temps discret.
Les paragraphes suivants sont dédiés à la présentation de leurs caractéristiques principales.
64
2. Si q(., .) est symétrique, ce qui est souvent le cas, alors les calculs sont plus simples
puisque α(xc , xt ) ne dépend que de la loi a posteriori :
π(xc )
α(xc , xt ) = min{ , 1}
π(xt )
Selon le choix de loi instrumentale, plusieurs versions de l’algorithme existent. On peut néan-
moins les classer en deux catégories :
ä soit la loi q(xc , xt ) est indépendante de xt , c’est la généralisation des approches d’accepta-
tion/rejet (Robert [60])
ä soit, le point xc est cherché dans un voisinage proche de xt : q(xc , xt ) = q(|xc − xt |)
Dans les deux cas, le choix implique une simplification de la probabilité d’acceptation et pré-
sente souvent de bonnes propriétés de convergence de l’algorithme.
Proposé par Geman et Geman en 1984 [22], l’algorithme de Gibbs est plus récent et à usage plus
spécifique que l’algorithme de Metropolis-Hastings. Il est particulièrement adapté aux struc-
tures multidimensionnelles de probabilités conditionnelles connues. Si x est le vecteur des in-
connues, dont la loi a posteriori est π(x), on suppose que x peut être partitionné en p sous
ensemblesnon vides (pas forcément composantes scalaires) : x = (x1 , x2 , ..., x p ) et on note :
ä x−i = x j , j 6= i le complémentaire de xi dans x.
ä π(xi |x−i ) = π(xi |x1 , . . . , xi−1 , xi+1 , . . . , x p ) la distribution conditionnelle de xi par rapport aux
composantes restantes.
L’échantillonnage par la méthode de Gibbs consiste à simuler une à une et de manière cyclique
chacune des composantes de x selon sa loi conditionnelle définie précédemment. Il peut être
résumé comme suit :
3. k ← k + 1 et reprendre l’étape 2
Remarque 1. Pour chaque composante i, la suite (xik )k est une chaîne de Markov.
Remarque 2. D’autre stratégies de visite existent comme par exemple le choix aléatoire à
l’itération k + 1 d’une composante à simuler, on l’appelle dans ce cas : échantillonneur de
Gibbs à balayage aléatoire.
65
Remarque 3. La convergence de l’algorithme est assurée si la loi cible vérifie la propriété de
positivité.
On constate qu’il n’est pas nécessaire de calculer la loi a posteriori de x mais qu’il suffit d’avoir
les lois conditionnelles sous-jacentes à sa partition, ce qui nous ramène à manipuler des lois
uni-variées. L’utilisation de la structure de la loi a posteriori pour la construction du noyau de
transition 3 , contrairement à l’algorithme de Metropolis-Hastings, présente l’avantage de géné-
rer des états avec un taux d’acceptation de 1.
Il peut parfois être compliqué d’évaluer une loi conditionnelle, on peut dans ce cas par exemple
utiliser l’algorithme de Metropolis-Hastings pour cette sous-étape ce qui conduit à une version
hybride de la méthode. Parmi les apports de la thèse, c’est justement le développement d’une
méthode hybride de Gibbs, appropriée aux structures récursives comme les séries chronolo-
giques et s’inscrivant dans le formalisme de filtrage de Kalman.
La section suivante sera dédiée à la mise en oeuvre de l’algorithme de Gibbs pour la caracté-
risation du bâtiment en temps discret. L’échantillonneur 4 de Kalman, bien qu’il soit d’aspect
théorique et général, sera présenté dans la partie traitant la simulation des températures (état
caché du bâtiment) afin de mieux cerner le besoin et le contexte de son application. L’étude du
problème dans le cas d’une discrétisation exacte, fortement non linéaire en les paramètres, peut
est être menée par l’adaptation de l’algorithme de Metropolis-Hastings.
Le choix de l’algorithme de Gibbs est motivé par ses grandes performances et par la propriété
de multilinéarité de la dynamique du bâtiment qui, dans le cadre gaussien adopté ici, simplifie
énormément les expressions des lois conditionnelles des inconnues.
4. Le terme « échantillonneur » est utilisé unanimement pour désigner le simulateur de Gibbs, nous l’utiliserons
également pour la simulation stochastique des séries temporelles
66
gaussiennes de paramètres facilement calculables. Comme nous allons le montrer au paragraphe
suivant, pour chaque itération de l’algorithme, la simulation se décline en trois étapes :
1. simulation des paramètres « constantes de temps » conditionnellement aux températures
(retenues de l’itération précédente) et aux paramètres « gains » selon la loi conditionnelle
π (Zt |XT , Zg ,Y ).
2. simulation des « gains » conditionnellement aux « constantes de temps » et aux tempéra-
tures suivant π (Zg |XT , Zt ,Y ).
3. Simulation des températures conditionnellement aux paramètres (générés aux étapes 1 et
2) selon la loi π (XT |Xθ ,Y ).
La formulation de chacune des lois conditionnelles utilisées précédemment est faite dans un
formalisme de filtre (et lisseur) de Kalman. L’intérêt de cette démarche est l’exploitation du
caractère récursif des sous-systèmes à simuler et des données dont on dispose. Le paragraphe
suivant est dédié à la présentation du filtre de Kalman dans le cas général et dans le contexte de
simulation (ou échantillonneur).
67
F IGURE 5.3: Schéma R3C2 pour le bâtiment
Comme nous le verrons, cela facilite grandement la mise en oeuvre de l’algorithme de Gibbs
dans le cadre gaussien.
Le bâtiment régulé en température peut être décrit par les équations suivantes :
1 1 1 1
Qres = Qch + AI = Cres .Ṫint + ( R f + Ri )Tint − Ri Ts − R f Text
1 1 1 1 (5.2)
Qs = Cs .Ṫs + ( Ro + Ri )Ts − Ri Tint − Ro Text
1
´
Qch = K E − Tint + τ (E − Tint )
où :
ä Qres : puissance de chauffage totale, elle est la somme des apports internes AI et de la part
puissance Qch fournie par le système de chauffage (convecteur, radiateur,...),
ä Qs : flux solaire,
ä E, Ts , Tint , Text : températures de consigne, de structure, intérieure et extérieure,
ä Cres , Cs , R1f , Ri
1 1
, Ro : capacités et conductances thermiques équivalentes du modèle du bâti-
ment,
ä K, τ : bande proportionnelle et constante de temps du régulateur Proportionnel Intégral de
chauffage.
Afin d’alléger les notations, on regroupe les paramètres inconnus dans un vecteur :
zres = C1res
zs = 1
C
z = 1s
τ τ
Z = z f = R1f et on rajoute une variable d’ état supplémentaire (équivalente à une tem-
zi = 1
R i
zo = 1
Ro
zK = K
pérature) : ˆ
1
d= (E − Tint )
τ
68
L’équation 5.2 de la boucle fermée devient :
zres Qres = Ṫint + zres (z f + zi )Tint − zres zi Ts − zres z f Text
zs Qs = Ṫs + zs (zo + zi )Ts − zs zi Tint − zs zo Text
d˙ = zτ (E − Tint )
Qch = zK (E − Tint + d)
E
ä observation (ou mesure) Qch .
L’écriture sous forme d’état déterministe de ce système peut alors se mettre sous la forme :
−zres (z f + zi + zK ) zres zi zres zK
Ẋ = zs zi −zs (zo + zi ) 0 X+
−zτ 0 0
zres z f zres 0 zres zK (5.3)
+ zs zo 0 zs 0 U
0 0 0 z
τ
Qch = zK −1 0 1 X+ 0 0 0 1 U
5.2.1.2 Discrétisation
Soit en discrétisant au pas de temps Te par la méthode d’Euler et en prenant en compte les bruits
d’état V et de mesure W :
(
X[k+1]−X[k]
= A (Zt , Zg ) X [k] + B (Zt , Zg )U [k] +V [k]
∆t (5.4)
Qch [k] = C (Zg ) X [k] + D (Zg )U [k] +W [k]
La méthode d’Euler (ou bloqueur d’ordre zéro) est une approximation qui se justifie dans notre
cas en comparant les valeurs propres de I + A∆t à celles de expA∆t .
Cette formulation du problème
estlinéaire en chacun des trois groupes d’inconnues :
Tint
ä Les températures X = Ts ,
d
5. En reportant Qres = Qch + AI
69
zres
ä Les « constantes de temps » Zt = zs ,
zτ
zf
zi
ä Les « gains » Z g = zo
zK
ce qui suggère de les simuler alternativement à l’aide de l’algorithme de Gibbs.
Le système 5.4 peut être réécrit de manière à permettre l’identification de la loi de Zt , en ex-
ploitant la linéarité. Soient les matrices :
− z f + zi + zK zi zK
M= zi − (zo + zi ) 0
−1 0 0
et
z f 1 0 zK
H = zo 0 1 0
0 0 0 1
et soit ΦT,Zg (X,U) la matrice diagonale, dont les coefficients sont les composantes du vecteur
MX + HU :
ΦT,Zg (X,U) = diag {MX + HU}
Alors l’itération décrite au 5.4 peut être réécrite :
X[k+1]−X[k]
∆t = ΦT,Zg (X [k] ,U [k]) .Zt +V [k]
v [k] étant supposé être une variable de loi gaussienne, ce qui conduit à la loi conditionnelle :
2
1 X[k + 1] − X[k]
log (π (X|Zt , Zg )) = cte − ∑ − ΦT,Zg (X [k] ,U [k]) .Zt (5.5)
2 k ∆t Q−1
70
2
1 X[k + 1] − X[k]
log (π (X, Zt , Zg )) = cte − ∑ − ΦT,Zg (X [k] ,U [k]) .Zt
2 k ∆t Q−1
+ log (π (Zt )) + log (π (Zg ))
La loi conditionnelle de Zt est identifiable en isolant les termes dépendant de Zt dans cette
expression :
1 1 1
log (π (Zt |X, Zg )) = cte − ZtT NZt + ZtT O + OT Zt + log (π (Zt ))
2 2 2
N
= ∑ΦT,Zg (X [k] ,U [k])T Q−1 ΦT,Zg (X [k] ,U [k])
k
avec :
O
= ∑ΦT,Zg (X [k] ,U [k])T Q−1 X[k+1]−X[k]
∆t
k
Ainsi, en adoptant une loi a priori gaussienne pour Zt :
1
log (π (Zt )) = cte − (Zt − µt )T Σt−1 (Zt − µt )
2
Nous obtenons la loi a posteriori pour Zt comme une gaussienne de covariance , ΣZt |X,Zg telle
que :
Σ−1 −1 T −1
Zt |X,Zg = Σt + ∑ΦT,Zg (X [k] ,U [k]) Q ΦT,Zg (X [k] ,U [k]) (5.6)
k
et d’espérance :
( )
X[k + 1] − X[k]
E [Zt |X, Zg ] = ΣZt |X,Zg µt + ∑ΦT,Zg (X [k] ,U [k])T Q−1 (5.7)
k ∆t
Nous reprenons le système 5.4 de la même manière que précédemment, en définissant cette
fois-ci :
AI
L = diag {Zt } Qs
E − Tint
et
Text − Tint Ts − Tint 0 E − Tint + d
ΦT,Zt (X,U) = diag {Zt } 0 Tint − Ts Text − Ts 0
0 0 0 0
Alors l’itération décrite au 5.4peut être réécrite :
71
X[k+1]−X[k]
∆t = ΦT,Zt (X [k] ,U [k]) .Zg + L [k] +V [k]
Lorsque l’on exprime la loi conditionnelle de Zg , le terme L [k] rentre dans la partie constante
et on retrouve les mêmes expressions de la covariance et l’espérance conditionnelles que 5.6 et
5.7.
Le conditionnement par rapport à la série de mesures complète est assuré par le lisseur, la mise
en oeuvre de ce dernier fournit les paramètres de la loi gaussienne π(Xk |Y1:N ) :
L’apport de la thèse à ce niveau consiste à générer de manière récursive une trajectoire X1:N S , si-
mulation de la loi conditionnelle π (X1:N |Y1:N ). Cette partie consacrée à l’échantillonneur achève
l’étape de simulation de X ∼ π(X|Zt , Zg ).
L’algorithme exploite la factorisation de la loi conditionnelle :
π (X1:N |Y1:N ) = π (X1 |X2:N ,Y1:N ) . . . π (XN−1 |XN ,Y1:N ) π (XN |Y1:N )
Remarque 5. Dans ce paragraphe ne figurent que les termes qui sont considérés comme aléa-
toires à cette étape de l’échantillonneur de Gibbs. Ni la commande, ni les paramètres Zt et Zg
n’apparaissent explicitement.
72
X0 ••• Xk Xk+1 ••• XN
Yk Yk+1 YN
5.2.3.1 Le filtre
73
La loi π(XN |Y1:N ) est obtenue en utilisant ces relations de proche en proche.
5.2.3.2 Le lisseur
Le lisseur de Kalman (voir Rauch et al [56] pour les équations) complète le filtre en condi-
tionnant par rapport à toutes les mesures, il procède de manière rétrograde à partir de XN . Les
espérance et variance conditionnelles de l’état à l’instant k :
µk|N = E {Xk |Y1:N }
Σk|N = cov Xk − µk|N
sont associées à la loi gaussienne π(Xk |Y1:N ). Pour établir les relations donnant les paramètres
de cette loi, on considère la loi conjointe de Xk et Xk+1 conditionnellement à toutes les mesures,
π(Xk , Xk+1 |Y1:N ) (en tenant compte des propriétés markoviennes de X1:N ) :
5.2.3.3 L’échantillonneur
S , on commence par la simulation conditionnelle de X , dont la loi est obte-
Pour simuler X1:N N
nue en fin de période du filtre. En effet, XN ∼ N µN|N , ΣN|N on peut donc prendre comme
échantillon :
74
S , X peut alors être échantillonné à partir de la loi π(X |X
de Xk+1 k k k+1 ,Y1:N ), dont les paramètres
peuvent être explicités à partir de 5.10 :
1h
π(Xk |Xk+1 ,Y1:N ) ∝ exp − (Xk+1 − AXk − BUk )T Q−1 (Xk+1 − AXk −
2
T i
−BUk ) + Xk − µk|k Σ−1
k|k Xk − µk|k
où Jk est le gain du lisseur de Kalman à l’instant k (voir annexe C). Par multiplication à gauche
par l’expression simplifiée de la variance on obtient l’espérance a posteriori :
Les équations 5.11 et 5.12 résument l’algorithme d’échantillonnage de Kalman, qui a permis de
S de l’état selon la loi conditionnelle π (X
générer une simulation X1:N 1:N |Zg , Zt ,Y1:N ,U1:N ). Ceci
achève la troisième étape d’une itération de l’algorithme de Gibbs.
Synthèse
Dans ce chapitre nous avons développé une version originale de l’algorithme de Gibbs, qui
permet de caractériser, par inversion bayésienne, le comportement thermique du bâtiment. Ce
−1 −1 T
7. A−1 + BCBT = A − AB C−1 + BT AB B A
75
dernier est modélisé sous forme d’un circuit R3C2 électrique équivalent asservi en température
par un régulateur Proportionnel Intégral.
L’approche que nous avons adoptée tire profit du caractère multilinéaire du système étudié et
fournit des simulations stochastiques des paramètres inconnus et de l’état. Pour générer les
simulations conditionnelles de l’état, nous avons mis au point un échantillonneur itératif basé
sur le formalisme de Kalman.
Cette démarche est généralisable dans le cadre d’un système stochastique gaussien, mis sous
forme de représentation d’état linéaire en l’état et multi-linéaire en les paramètres.
Par ailleurs, les travaux développés dans ce chapitre ont été implémentés sous forme d’une boite
à outils sous Matlab et ont fait l’objet d’une publication dans une conférence internationale
(Zayane et al [68]).
76
Chapitre 6
Ce chapitre est dédié à la présentation et l’analyse des résultats de l’algorithme de Gibbs appli-
qué au modèle global du bâtiment, agissant en boucle fermée. Ce modèle est constitué du circuit
électrique équivalent R3C2, pour l’enveloppe du bâtiment, couplé avec un régulateur Propor-
tionnel Intégral pour assurer le suivi de la température de consigne (température souhaitée).
La méthode a d’abord été appliquée à un jeu de données de synthèse, le choix des sollicitations
est fait avec un objectif de rapprochement avec le contexte réel. Les paramètres qui ont servi à
générer la courbe de charge sont connus et exploités dans un premier temps pour l’évaluation
des performances de la méthode. Dans un deuxième temps, nous nous intéresserons à l’étude
de la convergence de l’algorithme, et présenterons quelques outils graphiques et analytiques qui
permettent d’établir un diagnostic de convergence des distributions a posteriori.
L’approche d’inversion bayésienne est ensuite appliquée à des données obtenues par modéli-
sation fine d’un bâtiment constitué par cinq zones thermiques (pièces). La modélisation a été
réalisée sous le logiciel de description détaillée de bâtiments CLIM 2000, développé par EDF
et validé sur des données réelles.
Enfin des données réelles relevées sur site à partir d’un bâtiment du site de recherche des Renar-
dières d’EDF ont été exploitées pour déterminer les paramètres du bâtiment instrumenté. Ces
données mettent en évidence la nécessité de prendre en compte la limitation de la puissance
de chauffage qui se traduit par une transformation non linéaire du système global. L’effet de la
saturation sur ce dernier est étudié et nous proposons une manière de la modéliser dans un cas
simplifié où la régulation est assurée par une action proportionnelle.
77
zres ∗ coe fg ∗ coe ft
zres zs ∗ coe fg ∗ coe ft
coe fg = 1e4 zs
zτ ∗ coe ft
coe ft = 3600 zτ
z f coe fg
Qch ← Qch /coe fg zf ←
zi coe fg
Qs ← Qs /coe fg zi
zo coe fg
∆t ← ∆t/coe ft zo
zo coe fg
zk
zk coe fg
Tableau 6.1: Changement d’échelle des données et des variables
Ce changement d’échelle est rendu possible grace au fait que les paramètres interviennent par
produit. Il permet d’avoir une paramétrisation avec des conductances et capacités de même
ordre de grandeur sans que cela affecte les températures. Nous adopterons cette nouvelle para-
métrisation pour tous les traitements effectués dans la suite.
78
25
20
Températures (°C)
15
10
Tint
5 Ts
0 d
E
−5
0 200 400 600 800 1000 1200 1400 1600 1800 2000
Temps(s)
20
Qch
15
AI
Puissances (W)
10 Qs
−5
−10
0 200 400 600 800 1000 1200 1400 1600 1800 2000
Temps(s)
F IGURE 6.1: Données de synthèse : sollicitations, états et courbe de charge de chauffage asso-
ciés aux paramètres prédéfinis
Nous soulignons tout de même que ces passages en dessous de zéro sont ponctuels et pour la
plupart engendrés par le bruit de mesure, d’amplitude non négligeable.
Seules les puissances et les températures extérieure et de consigne serviront comme donnés
pour le problème d’inversion. Les résultats de l’algorithme seront présentés sous forme de si-
mulations stochastiques des sept paramètres scalaires et de l’état.
L’algorithme de Gibbs a été mis en oeuvre avec les mêmes variances des bruits d’état et de
mesure (que celles qui ont servi à générer les données de synthèse). Pour tester la méthode par
rapport à la connaissance a priori, nous avons d’abord opté pour le cas où elle n’est pas biaisée.
Ensuite, nous avons testé la méthode sur un jeu de paramètres avec des lois a priori gaussiennes
biaisées, ce qui traduit dans les cas de bâtiments réels l’incertitude par rapport à la connaissance
thermique ou sectorielle des valeurs des paramètres.
79
6.1.2.1 Connaissance a priori non biaisée des paramètres
La loi a priori associée à chaque paramètre est une gaussienne, de moyenne égale à la valeur
de référence (celle qui a servi à générer la courbe de charge) et d’écart type égal à la moitié de
celle-ci. re f
re f Z
ä Zg : µgpr = Zg , σgpr = diag( g2 )
re f
re f
ä Zt : µtpr = Zt , σtpr = diag( Zt
2 )
1
0.2
re f 2
ä µX0 = X0 et ΣX0 = diag( 0.1 )
0.2
Les simulations de la température intérieure sont illustrées sur la figure 6.2.
26
24
22
Tint(°C)
20
18
16
14
12
0 200 400 600 800 1000 1200 1400 1600 1800
4
3
(Tsint − Tint) (°C)
2
1
0
−1
−2
−3
0 200 400 600 800 1000 1200 1400 1600 1800
Temps (s)
Pour chaque itération de l’algorithme de Gibbs, une série temporelle Tints correspondant à T est
int
générée, elle est représentée par une couleur sur le premier graphe. Pour visualiser la variabilité
des températures simulées, l’écart entre chaque série et la température de référence (variable
d’état synthétisé) est illustré sur le deuxième graphe de la figure 6.2.
Les résultats de l’algorithme pour les paramètres sont les chaînes de Markov présentées sur la
figure 6.3.
Chaque graphe est associé aux simulations d’un des sept paramètres pour 3000 itérations de
l’algorithme de Gibbs. Les valeurs initiales des paramètres (initialisation de l’algorithme) sont
tirées aléatoirement à partir des lois gaussiennes a priori. L’analyse des simulations de la figure
6.3 montre qu’après une phase de démarrage (500 premières itérations), chaque chaîne atteint
un palier et varie autour de la valeur de référence du paramètre.
80
zres zk
4
0.24
0.22 3
0.2
0.18 2
0.16
1
0 500 1000 1500 2000 2500 3000 0 1000 2000 3000
zs zτ
0.15 0.4
0.1 0.3
0.05 0.2
0 500 1000 1500 2000 2500 3000 0 1000 2000 3000
zf zi zo
1 0.5
0.24
0.22 0.4
0.2 0.5
0.18 0.3
0.16
0 0.2
0 1000 2000 3000 0 1000 2000 3000 0 1000 2000 3000
Itération de GIBBS Itération de GIBBS Itération de GIBBS
F IGURE 6.3: Simulations des paramètres du bâtiment (noir) et du régulateur (bleu) pour 3000
itérations de l’algorithme de Gibbs, les valeurs de références sont tracées en rouge.
Afin de montrer l’apport de la méthode en terme de précision sur chaque paramètre par rapport
à la connaissance a priori, on peut par exemple comparer les fonctions de répartition relatives
aux lois a priori et a posteriori de chaque paramètre comme l’illustre la figure 6.4.
0.9
0.8
0.7
Probabilité cumulée
0.6
0.5 zpt
res
0.4
zpr
0.3 res
0.2
0.1
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
Données
F IGURE 6.4: Probabilités cumulées des lois a priori (magenta) et a posteriori (cyan) pour le
paramètre zres
81
Une chaîne de Markov est générée à partir de la loi a priori gaussienne et des fonctions de
répartitions gaussiennes sont ajustées aux chaînes correspondant aux distributions a priori et a
posteriori. La dispersion autour de la valeur de référence définit la vitesse de croissance de la
fonction de répartition. En effet, plus l’écart type de la loi gaussienne est petit, plus la croissance
de la fonction de répartition est rapide.
Les résultats de l’algorithme dans ce cas, où la connaissance a priori et les données mesurées
sont concordants, sont très satisfaisants. Dans la suite, on envisage le cas où la connaissance des
paramètres est moins précise dans le sens où les lois a priori sont biaisées.
Les espérances des loi a priori sont différentes des valeurs de référence, elles sont fixées à 0.5.
Les écart types sont choisis assez larges (de valeur 1) pour que la connaissance a priori ne soit
pas faussement trop informative, ce qui risque d’entraîner la non convergence de l’algorithme
ou une convergence trop lente.
ä Zg : µgpr = 0.5, σgpr = 1
ä Zt : µtpr = 0.5, σtpr = 1
Des exemples de simulations des capacités et des paramètres du régulateur sont illustrés sur la
figure 6.5.
zres zk
0.6
2.5
0.5
2
0.4
1.5
0.3
1
0.2
0.1 0.5
0 0
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
zs zτ
1
0.5 0.8
0.4
0.6
0.3
0.4
0.2
0.2
0.1
0 0
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
Iteration de Gibbs Iteration de Gibbs
F IGURE 6.5: Simulations des paramètres (bleu), espérances des lois a priori (rouge) et valeurs
de référence des paramètres en question (noir)
L’algorithme converge assez rapidement vers les valeurs de références pour ces quatre para-
mètres. Il reste néanmoins important d’étudier les propriétés de convergence de l’algorithme
82
plus en profondeur afin de s’assurer que les lois stationnaires sont bien explorées.
D’un point de vue pratique, il est aussi important d’estimer le nombre d’itérations nécessaires
pour l’algorithme pour décrire les lois stationnaires, pour réduire le temps de calcul. L’étude de
ces propriétés fait l’objet de la section suivante.
83
2.05
1.9
1.85
1.8
Burn in Régime stationnaire
1.75
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
Itération de l’algorithme MCMC
F IGURE 6.6: Exemple de simulation d’une chaîne de Markov par l’algorithme de Gibbs
Tracé des chaînes L’initialisation des paramètres est faite par un tirage aléatoire à partir de la
loi a priori, répéter cette expérience plusieurs fois et comparer les chaînes ainsi obtenues peut
être révélateur de plusieurs aspects :
ä La loi a posteriori peut avoir plusieurs modes (mélange gaussien par exemple), dans ce cas
on peut avoir des chaînes très différentes en fonction de l’initialisation. Le tracé des chaînes
représente, dans ce cas, un indicateur de l’homogénéité de la distribution a posteriori.
ä La longueur de la phase de burn in : en superposant les tracés des différentes chaînes, on
peut l’estimer plus précisément.
84
0.4 4
3
0.3
zres
zk
2
0.2
1
0.1 0
0 500 1000 1500 2000 0 1000 2000
0.4 0.5
0.3 0.4
zτ
zs
0.2 0.3
0.1 0.2
0 0.1
0 500 1000 1500 2000 0 1000 2000
0.4 0.8 1
0.3 0.6
zo
zf
zi
0.1 0.2
0 0 0
0 1000 2000 0 1000 2000 0 1000 2000
La figure 6.7 illustre des tracés de différentes chaînes des sept paramètres du modèle global du
bâtiment. La phase transitoire est par exemple plus visible pour les paramètres zres et zk . De
même, on peut déduire à partir de ces tracés que la variabilité de chaque chaîne est conservée
en changeant l’initialisation de l’algorithme, ce qui permet de pressentir une homogénéité des
lois a posteriori.
85
0.22
2.1
0.2 2.05
zres
zk
0.18 1.95
1.9
0.16
0.2
0.3
0.15
0.25
τ
zs
z
0.1
0.2
zo
zf
zi
0.15 0.4
0.2
0.1
0.2
0.1
F IGURE 6.8: Exemple de tracé sous forme de box-plot des simulations des sept paramètres
Ainsi ce tracé fournit une indication sur l’étendue de la chaîne (écart entre le minimum et
le maximum), une plage des valeurs les plus significatives (entre le premier et le troisième
quartile), ainsi qu’une valeur type (la médiane).
86
−4 −3
x 10 zres x 10 zk
1 3
2
0.5
1
0 0
0 50 100 150 200 250 300 350 400 450 500 0 100 200 300 400 500
−4 −3
x 10 zs x 10 zτ
4 1
2 0.5
0 0
0 50 100 150 200 250 300 350 400 450 500 0 100 200 300 400 500
−3 −3
x 10 zf x 10 zi zo
2 6 0.015
4 0.01
1
2 0.005
0 0 0
0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500
F IGURE 6.9: Variogrammes associés à un jeu de simulation des sept paramètres du modèle
global du bâtiment
portée » (distance nécessaire pour atteindre le palier) définit la distance au delà de laquelle les
valeurs prises par la chaîne ne sont plus liées. La portée est de ce fait un indicateur de la « durée »
de la période de démarrage (burn-in).
C’est une méthode proposée par Gelman et Rubin [21] en 1992. Elle consiste à analyser m
chaînes indépendantes, obtenues en lançant m fois l’algorithme de MCMC pour des initialisa-
tions différentes.
Cette heuristique permet de renseigner sur la vitesse de convergence (temps nécessaire pour
atteindre le régime stationnaire) et sur l’amélioration de la qualité des estimations, obtenues à
partir de la chaîne, en fonction du nombre d’itérations.
87
j=1:m
Supposons que la fonction à estimer est l’espérance, si chaque chaîne xi=1:N est de longueur
N = 2n, alors on ne prend en compte que les n derniers échantillons et on procède comme
illustré dans l’algorithme 6.1 :
1 m ¯j 2
B= ∑ x − x̄
m − 1 j=1
2n m
avec x¯j = ∑ x¯j
1 j 1
n ∑ xi et x̄ = m
i=n+1 j=1
2
2. On calcule la moyenne des m variances s j intra-séquences :
1 m j2
W= ∑s
m j=1
2 2n 2
avec sj = 1 j
∑ xi − x¯j
n−1
i=n+1
88
tenir compte de cette contrainte, nous présentons une méthode qui permet de générer un échan-
tillon d’un vecteur (de plusieurs composantes scalaires) suivant une loi gaussienne à support
positif.
Simulation d’une variable normale à support positif Plusieurs algorithmes existent (voir
Mazet et al [45] à titre d’exemple) pour la simulation d’une variable gaussienne scalaire tron-
quée :
!
(x − µ)2
f (x) = K exp − 1R+
2σ 2
1 si x ∈ I
où 1I est la fonction indicatrice : 1I (x) =
0 sinon
T
Simulation d’un vecteur à composantes gaussiennes tronquées Soit X = x1 · · · xl
un vecteur qui suit une loi gaussienne tronquée et que l’on cherche à simuler. Le principe de la
méthode consiste à appliquer un algorithme de Gibbs alternativement pour chacune des com-
posantes du vecteur à simuler.
Pour simuler X, il suffit pour cela d’évaluer les paramètres de la loi gaussienne de xi condition-
nellement aux autres composantes x−i .
89
3000 3000
zprior
res
zprior
k
2000 2000
zpost
res
zpost
k
1000 1000
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6
3000 2000
zprior
s 1500 zprior
τ
2000
zpost zpost
s τ
1000
1000
500
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.2 0.4 0.6 0.8
0 0 0
0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 0 0.5 1 1.5 2
90
2.4
0.22
zres
2.2
zk
0.2
2
0.18
0.2 0.35
0.15 0.3
zτ
zs
0.1
0.25
0.05
0.2
0.5
0.25
0.4 0.6
0.2
zo
zf
zi
0.3 0.4
0.15
0.2
0.1 0.2
0.1
0.5 0.25 0.2 0.5 0.25 0.2 0.5 0.25 0.2
F IGURE 6.11: Box-plot des lois a posteriori pour différents écart types des distributions a priori,
les écart types sont obtenus en multipliant l’espérance a priori respectivement par 0.5, 0.25 et
0.2
ment les variances des bruits d’état et de mesure. Les résultats de simulation sont illustrés sur
la figure 6.12.
Il ressort des tracés box-plot des lois a posteriori de la figure 6.12 une amélioration très nette
sur les paramètres lorsque les variances des bruits sont faibles, tant au niveau de l’espérance
(médiane autour de la valeur de référence) qu’au niveau de la précision (distance inter-quartiles
décroissante).
91
2.15
0.22 2.1
2.05
zres
zk
0.2
2
0.18 1.95
0.2 0.35
0.15 0.3
zτ
zs
0.25
0.1
0.2
0.1 0.5 1 0.1 0.5 1
0.5
0.25 0.6
0.4
0.2
zo
zf
zi
0.3 0.4
0.15
0.2
0.1
0.2
0.1
0.1 0.5 1 0.1 0.5 1 0.1 0.5 1
F IGURE 6.12: Quartiles des distributions a posteriori des paramètres en fonction des variances
des bruits d’état et de mesure, les matrices de variances de référence sont multipliées respecti-
vement par 0.1, 0.5 et 1. Les valeurs de référence sur les paramètres sont représentées par des
traits discontinus (noir).
Nous rappelons que le bâtiment modélisé est constitué de cinq zones thermiques soumises à la
même consigne de température (intermittence 19/20°C) et au même type de système de chauf-
fage, l’une des zones présente une saturation de la puissance de chauffage (courbe en magenta
sur la figure 6.13).
Le système équivalent considéré est supposé avoir une température intérieure, non disponible,
obtenue comme moyenne des cinq zones et une consommation égale à la somme des consom-
mations de toutes les zones.
92
E (°C) 20
19
0 10 20 30 40 50 60 70
8
Text (°C)
6
4
0 10 20 30 40 50 60 70
Qs (kW)
0.5
0
0 10 20 30 40 50 60 70
0.6
0.5
Qch (kW)
0.4
0.3
0.2
0.1
0 10 20 30 40 50 60 70
Temps (h)
F IGURE 6.13: Données CLIM2000 pour les 5 zones : température de consigne (E), température
extérieure (Text ), flux solaire (Qs ) et les courbes de charges correspondant au 5 zones.
ä la température de structure (Ts ) ne semble pas avoir une variabilité journalière, ce comporte-
ment n’est pas prévisible étant donné que toutes les sollicitations (températures de consigne
et extérieure, flux solaire) présentent varient à l’échelle de la journée.
Les grandeurs qui agissent directement sur la température de structure sont essentiellement les
données météo (Text et Qs ) et la température intérieure puisqu’elles agissent sur sa dérivée,
comme le montre bien la deuxième équation du système global 6.1 :
zres Qres = Ṫint + zres (z f + zi )Tint − zres zi Ts − zres z f Text
zs Qs = Ṫs + zs (zo + zi )Ts − zs zi Tint − zs zo Text
(6.1)
d˙ = zτ (E − Tint )
Qch = zK (E − Tint + d)
Or, comme on peut le lire à partir de la séquence des données d’entrée présentée sur la figure
6.13, la température extérieure ne varie que sur une plage de 4 °C et de manière assez lente
sur la durée d’observation (3 jours), ce qui atypique en pleine période de chauffe. De même,
l’amplitude des apports solaires est assez faible sur les 3 jours (cette amplitude n’est élevée
qu’en fin de période d’observation : troisième jours).
Quant à la température intérieure, ses variations sont définies par la température de consigne,
qui elle-même ne varie que sur 1°C. On peut donc conclure que, dans le cas où :
1. la météo est clémente (variations lentes et d’amplitude faible),
93
F IGURE 6.14: Simulations de l’état associé aux données CLIM2000, chaque couleur corres-
pond à une itération de l’algorithme de Gibbs, la température intérieure moyenne fournie par
CLIM2000 est reportée sur le premier graphe (noir).
2. le bâtiment est initialement « chargé » (un état initial avec absence de chauffage serait
sans doute plus intéressant),
3. le confort souhaité ne présente pas de variations considérables au cours de la journée,
la structure du bâtiment est dans un régime « pseudo-permanent » (pratiquement insensible aux
sollicitations).
Ce phénomène entraîne des difficultés considérables sur la démarche d’identification. L’analyse
des conséquences d’une telle configuration de bâtiment sur la précision des paramètres du mo-
dèle, avec un retour sur les résultats d’estimation du chapitre 4, fera l’objet de la suite de cette
section.
94
C’est en effet cette indétermination qui, combinée au fait que les paramètres interviennent par
produit sur les équations de la dynamique, explique la mauvaise qualité de la précision sur
l’estimation des paramètres, présentée au chapitre 4.
Il est ainsi évident que pour le cas des données CLIM2000, l’information sur le comportement
thermique du bâtiment ne se situe pas au niveau des valeurs des sept paramètres du modèle mais
à un niveau plus global. Et c’est justement pour de telles situations que la simulation présente
tout son intérêt. En effet, on peut facilement obtenir une approximation des distributions a
posteriori de plusieurs variables globales (gains statiques, résistance équivalente, produit de
variables) à partir des distributions des paramètres.
Sur la figure 6.15 présente le cas où la variable globale est le gain statique par rapport au flux
solaire.
40
GQs
35 a priori a posteriori
µ = 0.993428 || 0.986169 Gpr
Qs
σ = 0.124227 || 0.0500558
30
Histogramme
25
20
15
10
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Données
F IGURE 6.15: Histogrammes des distributions a priori (violet) et a posteriori (vert) du gain
statique par rapport au flux solaire
Le principe de l’approximation consiste, par exemple dans le cas où la variable globale est le
produit de deux paramètres, à multiplier les distributions des deux paramètres en question point
par point. On peut appliquer cette approximation pour les distributions a priori et a posteriori,
comme l’illustre le graphe 6.15.
95
6.5.1 Analyse des données et résultats
Les données, présentées sur la figure 6.16, sont échantillonnées au pas de temps de 10 minutes,
quand à la météo, elle est enregistrée sur le même site au niveau d’une cellule test (ETNA)
conçue pour l’expérimentation de nouveaux composants de bâtiment et la validation des logi-
ciels de modélisation. Le pas de mesure pour la météo est la minute et les données sont ensuite
sur échantillonnées au pas de 10 minutes.
LGEP1−10min−12/07
20
18
Consigne
T (°C)
16 Temperature interieure
14
12
10
500 1000 1500 2000 2500 3000 3500 4000
120
100
80
P (KW)
60
40 Puissance de chauffage
Apports internes
20
0
500 1000 1500 2000 2500 3000 3500 4000
Temps (10 minutes)
Le confort est défini avec une consigne de température intermittente c’est à dire en utilisant un
programmateur qui permet d’abaisser la température intérieure en période d’inoccupation tout
en maintenant la température de consigne désirée en période d’occupation.
96
Il n’y a pas d’amélioration significative des distributions a posteriori des paramètres par rap-
port aux distributions a priori, les histogrammes relatifs aux deux distributions, pour les sept
paramètres sont présentés sur la figure 6.17.
zres zK
1000
prior 1000
500 post
500
0
0 3.985 3.99 3.995 4 4.005 4.01 4.015
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
4
−8 x 10
x 10
zs zτ
600
600
400
400
200
200
0 0
0 0.5 1 1.5 2 2.5 3 3.5 1.1 1.105 1.11 1.115 1.12 1.125
−11 −3
x 10 x 10
zf zi zo
1500 1500 500
400
1000 1000
300
200
500 500
100
0 0 0
0 5000 10000 15000 0 0.5 1 1.5 2 2.5 0 1000 2000 3000 4000
4
x 10
F IGURE 6.17: Histogrammes des distributions a posteriori des paramètres du modèle du bâti-
ment LGEP1
On pourrait donc, à ce stade mettre en question la richesse des sollicitations, mais la simulation
de la température intérieure du bâtiment illustrée sur la figure 6.18 montre que, contrairement
au cas des données LGEP1 :
ä on ne retrouve pas la dynamique de cette composante de l’état (en particulier les descentes
de température sont en réalité beaucoup plus lentes) ;
ä le choix d’un réduit (valeur basse de la consigne) de 18°C n’est pas adapté à ce jeu de
données (biais sur le deuxième graphe).
Outre la richesse des sollicitations, d’autres causes peuvent être à l’origine du manque de perti-
nence de ces résultats, nous citerons notamment :
ä la limite de validité de l’approximation mono-zone pour le bâtiment LGEP1, ce qui remet-
trait en question le modèle R3C2 ;
ä la non prise en compte de la limitation de puissance, qui s traduirait par une saturation au
niveau du régulateur.
Pour confirmer l’existence de ces deux problèmes, nous avons enregistré la température inté-
rieure pour un bureau dans le bâtiment. Les mesures de température recueillies par plusieurs
sondes sont illustrées sur la figure 6.19.
97
s du bâtiment LGEP1, température
F IGURE 6.18: Simulations de la température intérieure Tint
moyenne Tint enregistrée (noir). L’écart entre chaque simulation et la moyenne est représenté
sur le deuxième graphe.
Une semaine de mesure est représentée sur cette figure. Le confort (valeur maximale de la
consigne) dans le bureau est de 22°C, nous constatons que :
ä les montées de température pour les quatre premiers jours de la semaine sont plus rapides
que la température moyenne du bâtiment, ceci est attribué au fait que les consignes ne sont
pas synchrones dans tous les bureaux. Le bâtiment ne peut donc être considéré comme une
zone thermique unique ;
ä le profil particulier du cinquième jour présente une montée de température plus lente que les
autres jours de la semaine, due à la présence de saturation.
Dans la suite de ce chapitre, nous nous intéressons à l’étude de la saturation, son effet sur
le comportement du bâtiment dans le cas où celui-ci est modélisé par un « R3C2+PI » et la
possibilité de l’intégrer dans l’algorithme de Gibbs.
98
Bureau LGEP1: semaine 4, 2009
26
25
sonde n° 1
24 sonde n° 2
pour différentes sondes
sonde n° 3
23 sonde n° 4
22
21
int
T
20
19
18
80 100 120 140 160 180 200 220 240
Temps (h)
F IGURE 6.19: Température intérieure d’un bureau appartenant au bâtiment LGEP1, pour une
période semblable (hiver)
chauffage est positive et majorée par une puissance Pmax , afin de prendre cette contrainte en
compte, deux approches peuvent être envisagées :
1. ne pas autoriser les saturations, ce qui consiste à obliger le système à être dans la partie
linéaire.
2. autoriser les saturations, ce qui fait apparaître des régimes de fonctionnement non li-
néaires du système en boucle fermée.
Or l’enjeu de la thèse n’est pas le choix du régulateur mais plutôt l’exploitation a posteriori des
données d’observations, le premier cas est donc à écarter et on se focalise dans la suite sur des
courbes de charges avec des zones de saturation. Soit u la puissance de chauffage requise par le
régulateur non saturé, la saturation consiste à remplacer u par 0 lorsqu’elle est négative et par
Pmax lorsqu’elle est supérieure à cette dernière, comme l’illustre la figure 6.20.
La commande saturée, effectivement injectée dans le système en boucle ouverte est :
Pmax si u(t) ≥ Pmax
sat(u(t)) = u(t) si 0 ≤ u(t) ≤ Pmax
0 sinon
99
sat(u(t))
Pmax
u(t)
0 Pmax
(E − Tint ) ce qui fait augmenter son intégrale en rampe. A la reprise de la phase linéaire de
fonctionnement (régulateur non saturé), le terme intégral est prédominant dans la correction et
risque d’entraîner de grandes variations de l’état du système.
Un test dans le cas du bâtiment, en prenant comme paramètres ceux qui ont servi à générer les
données de synthèse de la première partie du chapitre, est présenté sur la figure 6.21.
Sur cette figure, on remarque en particulier l’augmentation considérable de l’amplitude de l’ac-
tion intégrale lors des phases de saturation. Cette augmentation résulte en un dépassement de la
consigne qui peut atteindre 5 degrés et durer jusqu’à plus de la moitié du créneau (soit quelques
heures pour des cas réels).
Ces dépassements sont inacceptables parce qu’ils engendrent des situations d’inconfort d’une
part et parce qu’ils risquent d’aboutir à une instabilité de la boucle fermée. En pratique (dans
les systèmes réels), ce problème est traité par le ralentissement de la dynamique de l’action in-
tégrale en phase de saturation, cette action est appelée anti windup et fera l’objet du paragraphe
suivant.
Les modifications apportées au régulateur sont illustrées sur la figure 6.22 et résumées par le
système d’équations 6.2.
100
Action intégrale (°C)
d
20
dsat
−20
200 400 600 800 1000 1200 1400 1600 1800
−10
200 400 600 800 1000 1200 1400 1600 1800
Temps (s)
e u sat(u)
E + PI BO
- esat
+
+
- +
Ks
Tint
Qsat
ch = sat(Qch )
ż = e+Ks Qsat
ch −Qch = e+esat (6.2)
Qch = K e + 1 z
τi
Le principe de l’action anti windup est que esat atténue les variations de l’action intégrale en
régime saturé à travers son signe :
ä nul en régime linéaire : le fonctionnement du régulateur PI n’est pas affecté,
ä négatif quand Qch > Pmax : comme la puissance fournie Qsat ch est insuffisante pour réduire
l’écart positif entre E et Tint , e = E − Tint reste positif plus longtemps, ce qui fait augmenter
son intégrale z considérablement. esat , de signe négatif, permet de ralentir cette augmenta-
101
tion,
ä positif lorsque Qch < 0 : e est de signe négatif et le terme intégral diminue fortement. esat ,
de signe opposé à celui de e, permet de modérer la baisse de z.
Le poids attribué à l’action esat est pondéré par le gain de saturation Ks , qui représente un
paramètre d’ajustement par rapport à l’atténuation souhaitée de l’action intégrale. On notera
que le cas trivial où Ks est nul correspond à un système saturé sans action anti windup.
La figure 6.23 montre plusieurs profils de températures et de courbes de charge en fonction de
la valeur de Ks .
26 25
Ks=1
24 20
Ks=0.1
22 15 Ks=0.01
Tint (°C)
20 10
d (°C)
E
18 Ks=1 5
16 0
Ks=0.1
14 −5
Ks=0.01
12 −10
600 650 700 750 600 650 700 750
Temps (s) Temps (s)
11 8
10
6
9
4
Qch (W)
Ts (°C)
8 Ks=1
Ks=1 2
7 Ks=0.1
Ks=0.1
0
6 Ks=0.01
Ks=0.01
5 −2
600 650 700 750 600 650 700 750
Temps (s) Temps (s)
Il en ressort qu’un gain de saturation trop faible (0.01 par exemple) ne permet pas d’atténuer
suffisamment les dépassements de températures, et notamment celui de l’action intégrale, ce qui
fait durer la saturation plus longtemps. En revanche, un gain très élevé (1 par exemple) inverse la
dynamique de l’action intégrale (esat prépondérant par rapport à e) ce qui entraîne une chute de
la température intérieure et de la puissance de chauffage. Dans cette dernière configuration, le
système ne peut être considéré comme régulé par un PI avec saturation puisque la contribution
de ce dernier est complètement dominée par l’action anti windup.
Si le problème des dépassements engendrés par la saturation est résolu par la mise en oeuvre
d’une action anti windup, le ralentissement du régime transitoire du bâtiment reste une consé-
quence irrémédiable de la saturation. La figure 6.24 illustre l’écart entre les temps de réponse
du bâtiment avec ou sans saturation de la puissance de chauffage.
102
Action intégrale (°C)
6
d
dsat
4
Cette différence des temps de montée de la température intérieure s’explique par le fait qu’en
régime de saturation, le bâtiment agit en boucle ouverte. Il est soumis à une puissance de chauf-
fage constante égale soit à Pmax , si la consigne passe du réduit au confort, soit à 0 dans le cas
inverse.
103
ä l’équation de la dynamique décrit le sous système bâtiment, d’ordre deux (variables d’état
Tint et Ts ), soumis à trois sollicitations : Text , Qs et une entrée puissance : Qres = AI + Qch ,
ä l’équation d’observation décrit le sous système régulateur avec saturation, soit la dynamique
de la courbe de charge : Qch vue comme une fonction de la quantité (E − Tint ).
˙
Tint −zres (z f + zi ) zres zi 0 Tint
−zs (zo + zi ) 0 Ts +
Ts = zs zi
d −zτ 0 0 d
Text
zres z f zres 0 0
AI + Qch
+ zs zo 0 zs 0
Qs
0 0 0 zτ
E
Qch = S zK −1 0 1 X+ 0 0 0 1 U
et la boucle fermée qui décrit la phase linéaire, sans saturation et prend en compte le système
global : bâtiment et régulateur, avec la consigne comme entrée supplémentaire et la courbe de
charge comme observation (et non une sollicitation).
Le système mis sous cette forme a une matrice d’évolution dégénérée, ce qui pose des problèmes
d’instabilité. Pour s’affranchir de cette difficulté, deux possibilités peuvent être envisagées :
1. ne pas prendre en compte l’action intégrale, la régulation est assurée par une action pro-
portionnelle uniquement, on remarquera dans ce cas que la température de consigne n’in-
tervient qu’au niveau de l’équation d’observation,
2. amortir l’action intégrale (rajouter un paramètre pour la partie autonome de la dynamique
de ), ce qui peut avoir des conséquences considérables sur la dynamique du système global
et des justifications quant au choix de ce paramètre.
Synthèse
Dans ce chapitre, nous avons appliqué l’approche de simulation bayésienne pour identifier le
comportement thermique de bâtiment, modélisé par un système en boucle fermé « R3C2+PI »,
à plusieurs cas test.
L’étude des résultats correspondant au jeu de données de synthèse nous a permis de mettre
en évidence les performances de la méthode. Les chaînes simulées pour chaque paramètre du
modèle ont également servi à étudier la convergence de l’algorithme, et ce par l’application de
quelques critères graphiques.
Ensuite, nous avons appliqué l’algorithme d’inversion à des données obtenues par modélisation
fine d’un bâtiment constitué de cinq zones thermiques (pièces). L’analyse des simulations ob-
tenues a permis de déceler les causes de la mauvaise précision sur les paramètres obtenue au
chapitre 4, à savoir la dégénérescence du système due notamment à l’insuffisante variabilité de
la température de consigne et la clémence de la météo.
Par ailleurs, l’analyse des profils de température intérieure associés aux données réelles du
bâtiment de bureaux LGEP1 a montré les limites des approximations faites sur le modèle à
savoir la non prise en compte de la saturation et du caractère multi-zone du bâtiment. Nous
avons à cet effet étudié l’impact de la saturation sur le comportement thermique de bâtiment,
104
dans le cas où ce dernier est modélisé par un circuit électrique équivalent régulé en température
par un PI.
105
Conclusion Générale
Le secteur du bâtiment est le plus gros consommateur d’énergie en France avec des besoins
qui représentent 43% de l’énergie finale annuelle et qui correspondent à 25% des émissions de
CO2 au niveau national. Dans un contexte de préoccupation permanente d’économie d’énergie
et de conscience environnementale accrue, l’intérêt que présente le développement de stratégies
visant à minimiser la consommation d’énergie d’un bâtiment n’est plus à démontrer.
Que ces stratégies consistent à recommander l’amélioration de la structure du bâtiment, à sug-
gérer la modification des scénarios de gestion des systèmes de chauffage/climatisation ou à
préconiser certains comportements de l’usager, une démarche préalable de caractérisation du
comportement thermique de bâtiment s’avère inévitable.
La démarche consistant à identifier les paramètres d’un modèle de comportement thermique de
bâtiment, à partir de mesures recueillies sur site, n’est pas nouvelle et a été abordée :
ä soit sous l’angle du diagnostic de bâtiment : planification d’expériences le plus souvent
en régime libre, ou en régime régulé avec des hypothèses simplificatrices (régulateur de
dynamique très rapide par rapport au pas de mesure par exemple) ;
ä soit sous l’angle de l’élaboration de lois de commande pour le système de chauffage (ou de
climatisation) : tout ou partie des paramètres du modèle du bâtiment est fixée (identifiés au
préalable, connaissance thermique,...).
Le défi que nous nous sommes fixé pour ce travail est d’opter pour une démarche non-intrusive
pour l’usager, menée en mode d’occupation normale, c’est-à-dire en mode de régulation et qui
vise à identifier les modèles du bâtiment et du régulateur. Ce choix écarte la possibilité de
planifier une expérimentation du bâtiment et résulte naturellement en un ensemble de mesures
restreint incluant :
ä la température extérieure et le flux solaire global enregistrés à la station Météo France la
plus proche du bâtiment ;
ä la courbe de charge chauffage globale, relevée à partir d’un système de GTB (Gestion Tech-
nique du Bâtiment) ou un compteur intelligent et décomposée en deux parties : consomma-
tion de chauffage et apports internes ;
ä la température de consigne récupérée par GTB ou reconstruite par connaissance sectorielle.
Ces données constituent l’ensemble des entrées/sortie d’un système global « bâtiment + régula-
teur » composé :
1. d’un modèle R3C2 d’ordre 2 pour le bâtiment, proposé par EDF et obtenu par analogie
électrique/thermique. Ce modèle est de dynamique linéaire et il est non linéaire en les 5
paramètres qui le définissent (3 résistances et 2 capacités) ;
2. d’un modèle Proportionnel Intégral pour le régulateur.
107
Conscients de la difficulté qu’un cadre aussi contraint peut engendrer pour la démarche d’iden-
tification, notamment en l’absence de mesure de la température intérieure du bâtiment, nous
avons présenté les outils qui permettent d’apporter des éléments de réponse sur l’existence
d’une solution du problème et sur la qualité des estimations des paramètres identifiés.
Parmi ces outils, l’étude de l’identifiabilité nous a permis de vérifier que, pour le modèle glo-
bal choisi, l’unicité de la solution du problème d’identification est garantie d’un point de vue
théorique. Néanmoins, la précision que l’on peut espérer sur cette solution est fortement dé-
pendante de la « richesse » des sollicitations utilisées. Contrairement à l’identifiabilité, l’aspect
informationnel est plus qualitatif et les indicateurs que l’on peut utiliser, comme la matrice
d’information, nécessitent généralement des calculs complexes qui peuvent devenir rapidement
impraticables comme c’est le cas pour le système global que nous considérons.
C’est à ce titre que nous avons mené, dans une étape préliminaire à l’identification, une étude
de sensibilité globale et locale de la sortie du système (courbe de charge de chauffage) par
rapport à la variation des paramètres du modèle global afin de déterminer les directions ou
combinaisons des paramètres qui impactent le plus la sortie. Nous avons ensuite adopté une
démarche d’identification qui s’inscrit dans le cadre de l’estimation en boucle fermée, avec
comme inconnues le régulateur et le système en boucle ouverte et comme seule mesure(ou
sortie) la commande (grandeur qui assure la régulation).
Cette problématique d’identification par la commande pour un système en boucle fermée à
régulateur inconnu n’a, à notre connaissance, pas fait l’objet d’études antérieures.
Conformément à ce qui a été pronostiqué par l’étude de sensibilité, l’estimation des sept pa-
ramètres du modèle global à partir des données CLIM2000 (logiciel de simulation fine du
bâtiment) n’est pas précise, même si le modèle identifié permet de prédire la consommation
de chauffage de manière assez correcte. Dans l’objectif d’aller plus loin dans l’interprétation
de tels résultats et la prise en compte de connaissance a priori sur les paramètres, nous avons
proposé une deuxième approche qui consiste à solutionner le problème de l’identification sous
l’angle de la simulation stochastique, et ce dans un cadre bayésien.
Outre le fait qu’elle fournisse les distributions a posteriori des sept paramètres du modèle glo-
bal, la méthode que nous avons développée permet d’obtenir des simulations stochastiques des
composantes du vecteur d’état (températures intérieure et de structure notamment). Le principe
de cette approche consiste à utiliser les méthodes de MCMC et plus particulièrement l’échan-
tillonneur de Gibbs, pour générer les simulations des inconnues. L’échantillonneur de Gibbs
exploite la propriété mutilinéaire du système global par rapport aux trois groupes d’inconnues :
l’état, les paramètres « gains » (résistances et gain du régulateur) et les paramètres « constantes
de temps » (capacités et temps d’intégration du régulateur). Les lois conditionnelles associées
à chacun des groupes de paramètres sont d’expression explicite compte tenu du cadre gaussien
que nous avons adopté.
Pour simuler les séries temporelles (l’état du système), nous avons développé une méthode
qui exploite la théorie de lissage Kalman. En effet, le lisseur de Kalman est classiquement
utilisé pour calculer des espérances conditionnelles, notre apport consiste donc à adapter ce
formalisme au contexte de la simulation stochastique. La démarche que nous avons proposée
consiste à expliciter les lois conditionnelles et à proposer un algorithme qui permet de générer
des simulations de la série temporelle entière à partir de ces lois.
L’application de l’algorithme de simulation à l’identification du modèle global de comporte-
108
ment thermique de bâtiment est d’un apport considérable. A titre d’exemple, les simulations
de la température de structure pour les données CLIM2000 ont permis de déceler les causes
de l’indétermination, à savoir la dégénérescence du système due notamment à l’insuffisante va-
riabilité de la température de consigne et la clémence de la méteo. De même, nous avons pu
étudier le gain sur la précision apporté par la prise en compte d’informations supplémentaires
sur les paramètres.
Par ailleurs, l’analyse des profils de température intérieure associés aux données réelles du
bâtiment de bureaux LGEP1 a montré les limites des approximations faites sur le modèle à
savoir la non prise en compte de la saturation et du caractère multi-zone du bâtiment.
Pour pallier ces problèmes, nous proposons en guise de perspective à ce travail, la méthode de
simulation hamiltonienne qui permet d’intégrer ces contraintes, moyennant une linéarisation de
la saturation (par une fonction sigmoïde par exemple).
109
Approche hamiltonienne
La simulation hamiltonienne (voir à titre d’exemple Chen [11], Hanson [28] ou Neal [52]) est
une méthode de Monte Carlo hybride. L’idée de base de cette méthode est l’interprétation de la
loi a posteriori émanant de l’approche bayésienne comme l’exponentielle d’une énergie poten-
tielle d’un système physique. Cette énergie potentielle est augmentée d’une énergie cinétique,
l’ensemble formant ainsi un hamiltonien qui se conserve lors de l’évolution du système.
La nouvelle problématique consiste à simuler le hamiltonien, qui reste constant dans le cas
continu, et de choisir une méthode de discrétisation qui conserve les volumes pour que cette
propriété du hamiltonien soit conservée.
Principe
L’idée principale de cette méthode consiste à augmenter le vecteur x des paramètres, interprété
comme une « position », d’un vecteur y de même dimension vu comme une vitesse ou impul-
sion. La loi a posteriori à simuler est alors assimilée à une énergie potentielle U(x), elle est
étendue en hamiltonien par rajout d’une énergie cinétique relative aux impulsions :
1
H (x, y) = U (x) + |y|2
2
La chaîne de Markov (xt , yt ) est générée selon la probabilité stationnaire :
Mouvements dynamiques
Presque déterministes destinés à explorer les surfaces à H « constant » en respectant les équa-
tions du hamiltonien : (
dx ∂H
dt = ∂ y = y
dy
dt = − ∂∂Hx = −∇U (x)
Les simulations étant numériques, le choix d’un schéma de discrétisation adapté est très impor-
tant. Parmi les critères à respecter deux propriétés essentielles pour que l’état candidat respecte
l’équation détaillée :
111
ä la réversibilité temporelle : si la dynamique relie (x1 , y1 ) à (x2 , y2 ) alors elle relie (x2 , −y2 )
aussi à (x1 , −y1 )
ä la conservation des volumes dans l’espace des phases.
Ceci se conserve lorsqu’on discrétise par un schéma saute mouton (voir Hairer et al [27] pour
les schémas conservatifs) qui, appliqué au système hamiltonien précédent donne :
xi (t + τ) = xi (t) + τyi t + τ2
yi (t + τ) = yi t + τ2 − τ2 ∂∂Hx |x(t+τ)
(x̃t+1 , ỹt+1 ) si u ∈ U ([0, 1]) < exp (−∆H) ; ∆H = H (x̃t+1 , ỹt+1 ) − H (xt , yt )
(xt+1 , yt+1 ) =
(xt+1 , yt+1 ) sinon
Mouvements stochastiques
Ces mouvements concernent essentiellement la transition entre deux trajectoires hamiltonniennes
par génération du vecteur impulsion y. Ce choix peut être complètement aléatoire, en choisis-
sant un pas aléatoire sur une ou plusiers composante de y (voir [39] à titre d’exemple). mais
pour augmenter le taux de convergence de la méthode plusieurs variantes sont envisageables.
112
100 Leap Frog
Euler d ordre 4
90 Euler d ordre 1
80
Taux d acceptation en %
70
60
50
40
30
20
10
0
0 0.2 0.4 0.6 0.8 1
Pas de discrétisation, incrément de 0.001
Les inconnues du problème sont les impédances (résistances et capacités) ainsi que les séries
temporelles de température en quelques nœuds. Les équations de bilans des flux permettent
d’obtenir la partie énergie potentielle du hamiltonien.
113
Annexe A
Plan rez-de-chaussée Surface utile du niveau : 658,6 m², hauteur entre planchers : 3,18 m.
115
F IGURE A.2: Étage
Plan étage Surface utile du niveau : 630,7 m², hauteur entre planchers : 3,16 m.
Le système de chauffage fonctionne par intermittence du lundi au vendredi, la température de
confort (consigne haute) est autour de 21°C et le réduit (consigne basse) est de l’ordre de 18°C.
Les enregistrements sont faits au pas de 10 minutes à partir du mois de décembre 2007 et
comportent :
1. les températures intérieures de tous les locaux,
2. la consommation électrique totale du bâtiment,
3. les températures de consignes appliquées dans chacun des bureaux.
116
LGEP1−10min−12/07
20
18
Consigne
T (°C)
16 Temperature interieure
14
12
10
500 1000 1500 2000 2500 3000 3500 4000
120
100
80
P (KW)
60
40 Puissance de chauffage
Apports internes
20
0
500 1000 1500 2000 2500 3000 3500 4000
Temps (10 minutes)
LGEP1−10min−12/07
21.5
21
20.5
20
T (°C)
19.5
19
18.5
18
1750 1800 1850 1900 1950 2000
Temps (10 minutes)
117
LGEP1−10min−12/07
21.5 Consigne
Tint
21
20.5
T(°C)
20
19.5
19
18.5
60 Réelle
Filtrage morphologique
55
Seuillage
50
P (kW)
45
40
35
30
1030 1040 1050 1060 1070 1080 1090 1100 1110 1120 1130
Temps (10 min)
F IGURE A.6: Filtrage de la courbe de charge chauffage pour une durée de 16 heures
118
Courbe de charge chauffage LGEP1−10 min−12/07
120 Réelle
Seuillage
100 Filtrage morphologique
80
Qch (kW)
60
40
20
La fonction de transfert associée de ce modèle du premier ordre à trois entrées et une sortie est :
1 1 + aQ τs S
Tint (s) = Text (s) + Q(s) + ϕsh (s)
1 + τs GV (1 + τs) GV (1 + τs)
où s est la variable de Laplace.
Étant donné que la température intérieure est mesurée, on a choisi deux structures ARX de
modèles pour l’identification séparément de la régulation et du modèle du bâtiment comme
l’illustre la figure A.8.
La procédure a été mise en oeuvre séparément pour chacun des deux sous systèmes en considé-
rant comme jeux de données entrée/sortie :
1. Pour le régulateur : l’entrée est l’écart entre la température de consigne et la température
intérieure mesurée, la sortie est la courbe de charge chauffage
2. Pour le bâtiment : l’entrée est constituée de la météo (température extérieure et flux so-
laire) et de la somme entre la courbe de charge et les apports internes.
Un découpage de séquences de données servant à l’identification de chaque modèle à été adopté
en se basant sur des connaissances a priori fortes sur le comportement thermique du bâtiment à
savoir :
119
1. les montées de température (changement de la consigne du réduit à la température de
confort) traduisent essentiellement le fonctionnement du système de chauffage,
2. les descentes de la température intérieure (nuits) sont caractéristiques de la décharge des
parois du bâtiment, elles serviront par suite à identifier le modèle du bâtiment.
L’identification a été faite sous Matlab en prenant comme critère à minimiser l’erreur de pré-
diction, les résultats obtenus pour bâtiment, sont résumés dans le tableau A.1.
L’analyse des résultats d’identification ne permet pas de conclure à un modèle robuste du bâti-
ment et ce pour plusieurs raisons :
ä les résultats ne sont pas robustes par rapport au type de filtrage, mais aussi par rapport à la
séquence de données choisie (jours de la semaine),
ä les choix faits ne permettent pas de déterminer la surface sud équivalente puisque l’identifi-
cation est faite en l’absence d’ensoleillement,
ä une constante de temps du bâtiment de l’ordre 10 ou 15 heures n’est pas compatible avec la
connaissance thermique du bâtiment LGEP1,
ä les valeurs obtenues pour les déperditions statiques G du bâtiment ne permettent pas de
renseigner sur la structure de ce dernier. En effet les valeurs obtenues par identification pour
différentes configurations bornent les valeurs possibles pour n’importe quel bâtiment.
120
Base de données des paramètres de
bâtiment
121
Déperdition Rra (ren: 1 vol/h) Rp Rle Ro Ri Cs Cres Rhim Rhile Rheme Rhele
0.35W/m3.K 6,3E-04 INFINI 1,1E-03 1,3E-03 3,5E-05 452 169 103 140 173 567 1,91E-05 3,67E-04 3,60E-05 2,00E-04
0.35W/m3.K 6,3E-04 INFINI 7,2E-03 5,9E-04 7,5E-05 376 835 908 61 660 295 2,32E-05 2,29E-03 3,13E-05 1,25E-03
0.35W/m3.K 6,3E-04 INFINI 2,9E-03 7,0E-04 7,6E-05 356 571 803 61 611 365 2,35E-05 9,17E-04 3,25E-05 5,00E-04
0.35W/m3.K 6,3E-04 INFINI 1,1E-03 1,2E-03 7,7E-05 395 351 300 61 559 835 2,45E-05 3,67E-04 3,60E-05 2,00E-04
0.35W/m3.K 6,3E-04 INFINI 8,6E-04 2,0E-03 8,0E-05 393 264 831 60 728 518 2,50E-05 2,75E-04 3,83E-05 1,50E-04
0.35W/m3.K 6,3E-04 INFINI 7,2E-03 6,5E-04 2,1E-05 1 331 486 950 203 483 279 1,53E-05 2,29E-03 3,13E-05 1,25E-03
0.35W/m3.K 6,3E-04 INFINI 2,9E-03 7,5E-04 2,1E-05 1 347 443 990 199 989 673 1,54E-05 9,17E-04 3,25E-05 5,00E-04
0.35W/m3.K 6,3E-04 INFINI 8,6E-04 2,1E-03 2,3E-05 1 258 193 600 192 255 587 1,61E-05 2,75E-04 3,83E-05 1,50E-04
0.35W/m3.K 6,3E-04 INFINI 1,1E-03 1,3E-03 2,2E-05 1 289 636 140 196 095 606 1,58E-05 3,67E-04 3,60E-05 2,00E-04
0.36W/m3.K 6,3E-04 INFINI 7,2E-03 6,3E-04 3,3E-05 525 914 941 146 018 083 1,83E-05 2,29E-03 3,13E-05 1,25E-03
0.36W/m3.K 6,3E-04 INFINI 2,9E-03 7,4E-04 3,5E-05 405 710 734 145 607 463 1,85E-05 9,17E-04 3,25E-05 5,00E-04
0.35W/m3.K 6,3E-04 INFINI 8,6E-04 2,0E-03 3,6E-05 468 046 268 138 088 356 1,95E-05 2,75E-04 3,83E-05 1,50E-04
0.5W/m3.K 6,3E-04 INFINI 5,1E-03 3,9E-04 7,0E-05 420 037 178 64 944 358 2,32E-05 2,29E-03 3,13E-05 1,25E-03
0.5W/m3.K 6,3E-04 INFINI 2,0E-03 4,7E-04 7,2E-05 390 099 898 64 854 563 2,35E-05 9,17E-04 3,25E-05 5,00E-04
0.5W/m3.K 6,3E-04 INFINI 8,1E-04 8,2E-04 7,7E-05 327 345 490 64 575 122 2,45E-05 3,67E-04 3,60E-05 2,00E-04
0.5W/m3.K 6,3E-04 INFINI 6,1E-04 1,4E-03 7,8E-05 385 829 662 64 285 980 2,50E-05 2,75E-04 3,83E-05 1,50E-04
0.5W/m3.K 6,3E-04 INFINI 5,1E-03 4,5E-04 2,1E-05 1 310 683 360 205 330 141 1,53E-05 2,29E-03 3,13E-05 1,25E-03
0.5W/m3.K 6,3E-04 INFINI 2,0E-03 5,2E-04 2,1E-05 1 323 441 290 201 819 124 1,54E-05 9,17E-04 3,25E-05 5,00E-04
0.5W/m3.K 6,3E-04 INFINI 8,1E-04 8,8E-04 2,2E-05 1 267 440 640 199 029 940 1,58E-05 3,67E-04 3,60E-05 2,00E-04
0.5W/m3.K 6,3E-04 INFINI 6,1E-04 1,4E-03 2,2E-05 1 237 529 140 196 611 338 1,61E-05 2,75E-04 3,83E-05 1,50E-04
0.5W/m3.K 6,3E-04 INFINI 5,1E-03 4,3E-04 3,4E-05 452 244 457 147 113 781 1,83E-05 2,29E-03 3,13E-05 1,25E-03
0.5W/m3.K 6,3E-04 INFINI 2,0E-03 5,1E-04 3,4E-05 430 132 092 148 020 544 1,85E-05 9,17E-04 3,25E-05 5,00E-04
0.5W/m3.K 6,3E-04 INFINI 8,1E-04 8,6E-04 3,5E-05 389 161 124 149 345 312 1,91E-05 3,67E-04 3,60E-05 2,00E-04
0.5W/m3.K 6,3E-04 INFINI 6,1E-04 1,4E-03 3,5E-05 461 729 189 142 735 439 1,95E-05 2,75E-04 3,83E-05 1,50E-04
0.8W/m3.K 6,3E-04 INFINI 5,7E-04 4,2E-04 7,1E-05 355 378 515 71 076 660 2,45E-05 3,67E-04 3,60E-05 2,00E-04
0.8W/m3.K 6,3E-04 INFINI 4,3E-04 4,5E-04 7,2E-05 328 216 424 73 339 880 2,50E-05 2,75E-04 3,83E-05 1,50E-04
0.8W/m3.K 6,3E-04 INFINI 3,6E-03 2,7E-04 2,1E-05 1 265 046 100 207 847 526 1,53E-05 2,29E-03 3,13E-05 1,25E-03
0.8W/m3.K 6,3E-04 INFINI 1,4E-03 3,1E-04 2,1E-05 1 279 673 480 205 019 016 1,54E-05 9,17E-04 3,25E-05 5,00E-04
0.8W/m3.K 6,3E-04 INFINI 5,7E-04 4,8E-04 2,2E-05 1 238 408 330 203 763 363 1,58E-05 3,67E-04 3,60E-05 2,00E-04
0.8W/m3.K 6,3E-04 INFINI 4,3E-04 6,8E-04 2,2E-05 1 215 611 140 203 314 356 1,61E-05 2,75E-04 3,83E-05 1,50E-04
0.8W/m3.K 6,3E-04 INFINI 3,6E-03 2,6E-04 3,0E-05 411 854 025 180 950 424 1,83E-05 2,29E-03 3,13E-05 1,25E-03
0.8W/m3.K 6,3E-04 INFINI 1,4E-03 3,0E-04 3,0E-05 398 467 983 180 492 422 1,85E-05 9,17E-04 3,25E-05 5,00E-04
0.8W/m3.K 6,3E-04 INFINI 5,7E-04 4,7E-04 3,2E-05 345 150 788 178 070 319 1,91E-05 3,67E-04 3,60E-05 2,00E-04
0.8W/m3.K 6,3E-04 INFINI 4,3E-04 6,6E-04 3,3E-05 321 250 416 179 765 397 1,95E-05 2,75E-04 3,83E-05 1,50E-04
Annexe B
1 N−1
log [ fθ (Y )] = cste − ∑ kyk+1 − ϕθ (y0:k )k2Σ−1
2 k=0 ε
123
et par suite sa hessienne est :
1 N−1 n T 2T
o
∇2θ log [ fθ (Y )] = − ∑ ∇θ ϕθ (y0:k ) Σ−1
ε ∇θ ϕ θ (y0:k ) − ∇θ ϕ θ (y0:k ) Σ −1
ε (yk+1 − ϕ θ (y0:k ))
2 k=0
T
or Eθ [yk+1 ] = ϕθ (y0:k ) et ∇2θ ϕθ (y0:k )est une constante à l’instant k + 1 donc :
h T i
2 −1
Eθ ∇θ ϕθ (y0:k ) Σε (yk+1 − ϕθ (y0:k )) = 0 (B.3)
1 N−1 T
Eθ ∇θ ϕθ (y0:k ) Σ−1
I (θ ) = ∑ ε ∇θ ϕθ (y0:k )
2 k=0
En effet, en notant
∂ (yk − ŷk ) ∂ ek (θ )
ψ(k, θ ) = − =−
∂θ ∂θ
N
∂ J(θ )
le gradient de J est ∂θ = − N2 ∑ ek (θ )ψ T (k, θ )
k=1
N N
∂ 2 J(θ ) 2 ∂ 2 ek (θ )
La hessienne est : ∂θ2
= N ∑ ψ(k, θ )ψ T (k, θ ) + N2 ∑ ek (θ ) ∂θ2
k=1 k=1
Dans le cas où la structure du modèle est fonction d’une paramétrisation fixée, la hessienne
risque d’être singulière. Ceci peut entraîner des cas de non convergence du fait que la plupart
des algorithmes de minimisation utilisent l’inverse de la hessienne. Plusieurs régularisations
existent dont :
124
ä Gauss Newton : plus spécifique à la minimisation de somme de carrés, il présente l’avantage
de ne pas nécessiter des dérivées secondes. La dérivée d’ordre 2 est négligée dans le calcul
de la hessienne :
" #T " #
∂ 2 J(θ̂ i ) 2 N ∂ Jr (θ̂ i) ∂ Jr (θ̂ i)
≈ ∑ ψ(k, θ )ψ T (k, θ ) = 2
∂θ2 N k=1 ∂θ ∂θ
e1 (θ̂ i )
..
où Jr est la matrice jacobienne du résidu : r = calculée en θ̂ i .
.
eN (θ̂ i )
ä Levenberg Marquardt : c’est une version relaxée de Gauss Newton :
" #T " #
∂ 2 J(θ̂ i ) ∂ Jr (θ̂ i ) ∂ Jr (θ̂ i )
≈2 +λI
∂θ2 ∂θ ∂θ
125
Annexe C
Filtrage de Kalman
Le filtre de Kalman répond au problème d’estimation en s’appuyant sur la notion d’état. l’intérêt
de la représentation d’état est qu’elle permet de caractériser des systèmes invariants et non
invariants et donc de modéliser des signaux aléatoires stationnaires ou pas.
Comme on se place dans un cadre de données discrètes (courbe de charge et données météo),
nous présenterons le filtre dans le cas de système discret.
127
ä E VkVl = Qk δkl , Q est la matrice de covariance de V .
T
ä E WkWl = Rk δkl , R est la matrice de covariance de W .
T
C.2 Le filtre
Le principe du filtrage linéaire récursif de Kalman est :
1. Prédire l’état à l’instant k + 1 à partir de son estimation à l’instant k. Soit en notant µk|k
l’estimation de X à l’instant k, et µk+1|k sa prédiction à partir de l’estimation à l’instant k,
on a :
µk+1|k = Ak µk + BkUk
Cette prédiction n’est pas bruitée car il s’agit d’une projection sur l’espace formé par les
observations et l’état jusqu’à l’instant k et que par définition le bruit d’état est blanc (les
valeurs qu’il prend aux différents instants ne sont pas corrélées, il est par suite non corrélé
avec l’état aux instants précédents) et non corrélé avec le bruit de mesure (donc avec les
observations). Soit Σk+1|k la variance de l’erreur de prédiction :
T
E Xk+1 − µk+1|k Xk+1 − µk+1|k
T
et Σk|k la variance de l’erreur d’estimation (E Xk − µk|k Xk − µk|k ), on peut facile-
ment montrer que :
Σk+1|k = Ak Σk|k ATk + Qk
2. Estimer après l’« arrivée » de la mesure Yk+1 : cette phase permet de corriger la prédiction
en ajoutant un terme correctif (écart entre la mesure et son estimation) pondéré par un
gain K (de Kalman) tel que :
µk+1|k+1 = µk+1|k + Kk+1 Yk+1 −Ck+1 µk+1|k
T −1
Ck+1 Σk+1|kCTk+1 + Rk+1
Kk+1 = Σk+1|kCk+1
Σk+1|k+1 = {I − Kk+1Ck+1 } Σk+1|k
C.3 Le lisseur
Le principe du lisseur de Kalman est de mettre en oeuvre le filtre dans un premier temps et en-
suite de faire le sens inverse en commençant par la dernière estimation et en allant récursivement
jusqu’à la première. En conservant les notations présentées précédemment et en posant :
ä µk+1|k la prédiction dans le sens inverse,
ä Σk+1|k la variance de l’erreur de prédiction correspondante,
128
ä Σk|N la variance de l’erreur d’estimation du lisseur et Jk le gain du “retour”.
Les équations à mettre en oeuvre à l’itération k du sens inverse sont :
129
Bibliographie
[2] A. A NDROUTSOPOULOS, J.J. B LOEM, H.A.L. van D IJK et P.H. BAKER : Comparison
of user performance when applying system identification for assessment of the energy
performance of building components. Building and environment, 43:189–196, 2008. 2.2.1
[3] P. R. A RMSTRONG : Model identification with application to building control and fault
detection. Thèse de doctorat, Massachusetts Institute of Technology, 2004. 2.2.2.2
[4] A. A. El A ZHER : Identification des systèmes à deux échelles de temps et application au
chauffage optimal des bâtiments. Thèse de doctorat, Ecole des Mines de Paris, 1992. 1.3,
2.2.2.2
[5] S. BARROIS et J.M J ICQUEL : Fdd usages : Outil de pré-diagnostic automatique à partir
de la courbe de charge annuelle d’un site tertiaire. Rapport technique, EDF, 2003. 1.2.2.1
[6] A. S. BAZANELLA, M. G EVERS et L. M ISKOVIC : Closed-loop identification of mimo
systems : a new look at identifiability and experiment design. In Proceedings of European
Control Conference, 2007. 3.2.1
[7] I. B RAEMS, L. JAULIN, M. K IEFFER et E. WALTER : Guaranteed numerical alternatives
to structural identifiability testing. Prodeedings of the 40th IEEE Conference on Decision
and Control, 4:3122–3127, 2001. 2
[8] S. P. B ROOKS et G. O. ROBERTS : Assessing convergence of markov chain monte carlo
algorithms. Statistics and Computing, 8:319–335, 1997. 6.2
[9] M.J. C HAPPELL : Structural identifiability and indistinguishability of certain two-
compartment models incorporating nonlinear efflux from the peripheral compartment.
Mathematical Biosciences, 125:61–81, 1995. 3.2, 1
[10] M.J. C HAPPELL et R.N. G UNN : A procedure for generating locally identifiable reparame-
trisations of unidentifiable non-lineai systems by the similarity transformation approach.
Mathematical Biosciences, 148:21–41, 1998. 1
[11] L. C HEN, Z. Q IN et J.S. L IU : Exploring hybrid monte carlo in bayesian computation.
ISBA and Eurostat, 2001. 6.6.2.1
[12] D. C OVALET et N. P ILLA : Modélisation dynamique sur ca-sis version 3.0 du bâtiment de
bureaux lgep1. Rapport technique, EDF, 2005. A.1
131
[13] S. DAUTIN : Réduction de modèles thermiques de bâtiments : amélioration des techniques
par modélisation des sollicitations météorologiques. Thèse de doctorat, Université de
POITIERS, 1997. 2.1.2.2, 2.1.2.3, 2.2.1
[14] T.J. D EBUS, P.E. D UPONT et R.D. H OWE : Distinguishability and identifiability testing
of contact state models. Advanced Robotics, 19:545–566, 2005. 3.2
[15] F. D ÉQUÉ : Techniques numériques de réduction de modèles linéaires : état de l’art.
Rapport technique, EDF, 1994. 2.1.2.3
[16] W. F ELLER : An Introduction to Probability Theory and Its Applications. John Wiley,
1968. 5.1.2
[17] U. F ORSSELL : Properties and Usage of Closed-loop Identification Methods. Thèse de
doctorat, Linkoping University, 1997. 3.3.2
[18] G. F RAISSE : La Régulation Thermique Des Bâtiments Tertiaires : Application De La
Logique Floue A La Régulation Centrale Du Chauffage En Régime Intermittent. Thèse de
doctorat, Institut National Des Sciences Appliquées De Lyon, 1997. 1.3, 2.2.2.2
[19] Y. G AO, J.J. ROUX, C. T EODOSIU et L.H. Z HAO : Reduced linear state model of hollow
bloks walls, validation using hot box measurements. Energy and Buildings, 36:1107–
1115, 2004. 2.1.2.3
[20] A. E. G ELFAND et A. F. M. S MITH : Sampling-based approaches to calculating marginal
densities. Journal of the American Statistical Association, 85:398–409, 1990. 5.1.2.1
[21] A. G EMAN et D. B. RUBIN : Inferences from iterative simulating using multiple se-
quences. Statistical Science, 7:457–511, 1992. 6.2.2.1
[22] S. G EMAN et D. G EMAN : Stochastic relaxation, gibbs distributions, and the bayesian
restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence,
6:721–741, 1984. 5.1.2.3
[23] M. G EVERS, A.S. BAZANELLA, X. B OMBOIS et L. M ISKOVIC : Identification and the
information matrix : how to get just sufficiently rich ? IEEE Transactions on Automatic
Control, 54:2828–2840, 2009. 3.2.1
[24] M. G EVERS, L. M ISKOVIC, D. B ONVIN et A. K ARIMI : Identification of multi-input
systems : varaiance analysis and input design issues. Automatica, 42:559–572, 2006. 1, 2
[25] H. G HARBI, M. G ONORD et A. M ARTI : Reconnaissance de la signature thermique d’un
logement. signa/ther. Rapport technique, EDF, 1989. 2.1.2.1
[26] W.R. G ILKS, S. R ICHARDSON et D.J. S PIEGELHALTER : Markov Chain Monte Carlo In
Practice. Chapman & Hall, 1995. 5.1.2.1, 2
[27] E. H AIRER, C. L UBICH et G. WANNER : Geometric Numerical Integration. Springer,
2004. 6.6.2.1
[28] K.M. H ANSON : Markov chain monte carlo posterior sampling with the hamiltonian
method. Medical Imaging : Image Processing, 4322:456–467, 2001. 6.6.2.1
132
[29] W. K. H ASTINGS : Monte carlo sampling methods using markov chains and their appli-
cations. Biometrica, 57:97–109, 1970. 5.1.2.2
[30] P. M. J. Van Den H OF et R. J. P. S CHRAMA : Identification and control : closed loop
issues. Automatica, 31:1751–1770, 1995. 3.3.2
[31] S. JAZOULI et S. BARROIS : Calage d’un modèle thermique de bâtiment "r3c2" à partir
de la courbe de charge et de la météo. Rapport technique, EDF, 2002. 2.3
[32] A.H. JAZWINSKI : Stochastic Processes And Filtering Theory. Academic Press, 1970.
5.2.3.1
[33] M.J. J IMENEZ et H.M ADSEN : Models for describing the thermal characteristics of buil-
ding components. Building and Environment, 43:152–162, 2008. 2.2.1
[34] J. K AIPIO et E. S OMERSALO : Statistical and Computational Inverse Problems. Springer,
2004. 5.1.1
[35] R. E. K ALMAN : A new approach to linear filtering and prediction problems. Transactions
of the ASME–Journal of Basic Engineering, 82:35–45, 1960. 5.2.3.1
[36] S. L AGRANGE, N. D ELANOUE et L. JAULIN : Injectivity analysis using interval analysis
application to structural identifiability. Automatica, 44:2959–2962, 2008. 2
[37] I. D. L ANDAU et A. B ESANÇON -VODA : Identification des systèmes. Hermes Science,
2001. 1, 2, B.2
[38] L. L ARET : Contribution au développement de modèles mathématiques du comportement
thermique transitoire des structures d’habitation. Thèse de doctorat, Université de LIEGE,
1980. 2.1.2.3
[39] B. L EIMKUHLER et S. R EICH : Simulating Hamiltonian Dynamics. Cambridge University
Press, 2005. 6.6.2.1
[40] L. L JUNG : Asymptotic variance expressions for identified black-box transfer function
models. IEEE Transactions on Automatic Control, 30:834–844, 1985. 2
[41] L. L JUNG : System Identification - Theory For the User. Prentice Hall ; 2 edition, 1999.
3.2.1
[42] S. A. M ARSHALL : An approximate method for reducing the order of a linear system.
Control, 10:642–643, 1966. 2.1.2.3
[43] A. M ARTI et J.P. R IGNAC : Méthode de résolution du système différentiel utilise par le
simulateur thermique de stratege 2.1. Rapport technique, EDF, 1991. 2.1.2.3, 2.3
[44] G. M ATHERON : The theory of regionalized variables and its applications. Les cahiers du
CMM de Fontainebleau, Ecoles des Mines de Paris, 1971. 6.2.1
[45] V. M AZET, D. B RIE et J. I DIER : Simuler un distribution normale à support positif à partir
de plusieurs lois candidates. In GRETSI, 2005. 6.2.3
[46] E. M AZRIA : Le guide de l’énergie solaire passive. Parenthèses, 1981. (document), 1.4
133
[47] T. M C K ELVEY : Identification of State-Space Models from Time and Frequency Data.
Thèse de doctorat, Linkoping University, 1995. 4.2.3
[48] C. M ENEZO : Contribution à la modélisation du comportement thermique des bâtiments
par couplage de modèles réduits. Thèse de doctorat, INSA de Lyon, 1999. 2.1.2.3
[49] N. M ETROPOLIS et A. W. ROSENLUTH : Equation of state calculations by fast computing
machines. Journal of Chemical Physics, 21:1087–1092, 1953. 5.1.2.1, 5.1.2.2
[50] B. C. M OORE : Principal component analysis in linear systems : controllability, observa-
bility, and model reduction. IEEE Transactions on Automatic Control, 26(1):17–32, 1981.
2.1.2.3
[51] Direction de la Climatologie M ÉTÉO F RANCE : Fiche méthode degrés jours. Dis-
ponible sur internet à l’adresse http ://climatheque.meteo.fr/Docs/DJC-methode.pdf, 03
2005. 2.1.1.1
[52] R.M. N EAL : Probabilistic inference using markov chain monte carlo methods. Rapport
technique, University of Toronto, 1993. 6.6.2.1
[53] B. N INNESS et H. H JALMARSSON : Variance error quantifications that are exact for finite
model order. IEEE Transactions on Automatic Control, 49:1275–1291, 2004. 2
[54] F. O LLIVER et A. S EDOGLAVIC : Algorithmes efficaces pour tester l’identifiabilité lo-
cale. In Actes de la Conférence Iternationale Francophone d’Automatique, pages 811–
816, 2002. 2
[55] D. P ETIT : Réduction de modèles de connaissance et identification de modèles d’ordre
réduit. Thèse de doctorat, Université de Provence (Aix-Marseille I), 1991. 2.1.2.3
[56] H. E. R AUCH, F. T UNG et C. T. S TRIEBEL : Maximum likelihood estimates of linear
dynamic systems. AIAA J., 3:1445–1450, 1965. 5.2.3.2
[57] J. C. R IBOT, A. H. ROSENFELD, W. L UHRSEN et F. F OUQUET : What Works : Documen-
ting the Results of Energy Conservation in Buildings, chapitre Monitored Low-Energy
Houses in North America and Europe : A Compilation and Economic Analysis. American
Solar Energy Society, 1983. 2.1.1.1
[58] V. R ICHALET : Caracterisation energetique des bâtiments sur site. Identification de mo-
deles dynamiques. Methodes de signature Energetique. Thèse de doctorat, INPG, 1991.
1.3, 2.1.1, 2.2.2.1
[59] C. ROBERT : L’analyse statistique Bayesienne. Economica, 1992. 5.1.1
[60] C. ROBERT : Méthodes de Monte Carlo par chaînes de Markov. Economica, 1996. 5.1.2.1,
5.1.2.2, 5.1.2.2
[61] C. P. ROBERT : Discretization and MCMC Convergence Assessment. Springer, 1998. 6.2
[62] J.J. ROUX : Proposition de modèles simplifiés pour l’étude du comportement thermique
des bâtiments. Thèse de doctorat, INSA Lyon, 1984. 2.1.2.3
[63] T. S ODERSTROM et P. S TOICA : System Identification. Prentice Hall, 1989. 3.1.2, 4.2
134
[64] E. WALTER, Y. L ECOURTIER et J. H APPEL : On the structural output distinguishability
of parametric models, and its relations with structural identifiability. IEEE Transactions
on Automatic Control, 29-1:56–57, 1984. 3.2
[65] E. WALTER et L. P RONZATO : Identification de modèles paramétriques à partir de don-
nées expérimentales. Masson, 1994. 3.1.2, 3.2
[66] E. WALTER et L. P RONZATO : On the identifiability and distinguishability of nonlinear
parametric models. Mathematics and Computers in Simulation, 42:125–134, 1996. 1
[67] B. Y U et P. M YKLAND : Looking at markov samplers through cusum path plots : A simple
diagnostic idea. Statistics and Computing, 8(3):275–286, 1998. 6.2.2.2
[68] C. Z AYANE, C. L AJAUNIE, L. P RALY, A. G IRARD et J. M. J ICQUEL : Joint state and
parameter estimation : a bayesian approach. IEEE ICMSC, 1:137–142, 2010. 5.2.3.3
135
Identification d’un modèle de comportement thermique de bâtiment à partir de
sa courbe de charge
Résumé : Dans un contexte de préoccupation accrue d’économie d’énergie, l’intérêt que présente
le développement de stratégies visant à minimiser la consommation d’un bâtiment n’est plus à dé-
montrer. Que ces stratégies consistent à recommander l’isolation des parois, à améliorer la gestion
du chauffage ou à préconiser certains comportements de l’usager, une démarche préalable d’identifi-
cation du comportement thermique de bâtiment s’avère inévitable.
Contrairement aux études existantes, la démarche menée ici ne nécessite pas d’instrumentation du
bâtiment. De même, nous considérons des bâtiments en occupation normale, en présence de régula-
teur de chauffage : inconnue supplémentaire du problème. Ainsi, nous identifions un système global
du bâtiment muni de son régulateur à partir de :
1. données de la station Météo France la plus proche ;
2. la température de consigne reconstruite par connaissance sectorielle ;
3. la consommation de chauffage obtenue par système de Gestion Technique du Bâtiment ou par
compteur intelligent ;
4. autres apports calorifiques (éclairage, présence de personnes...) estimés par connaissance sec-
torielle et thermique.
L’identification est d’abord faite par estimation des paramètres (7) définissant le modèle global, en
minimisant l’erreur de prédiction à un pas. Ensuite nous avons adopté une démarche d’inversion
bayésienne, dont le résultat est une simulation des distributions a posteriori des paramètres et de la
température intérieure du bâtiment.
L’analyse des simulations stochastiques obtenues vise à étudier l’apport de connaissances supplé-
mentaires du problème (valeurs typiques des paramètres) et à démontrer les limites des hypothèses
de modélisation dans certains cas.
Mots clés : Identification, inversion bayésienne, échantillonneur de Gibbs, simulation stochastique
Identification of a building thermal behavior model from its load curve
Abstract: In a context of permanent concern of energy saving, the importance of developing strate-
gies to minimize the energy consumption of a building is not to be any more demonstrated. Whether
these strategies consist in recommending building insulation, suggesting the modification of the man-
agement of heating/air conditioning systems or advising the change of certain occupants’ habits, a
preliminary step of characterizing the thermal behavior of the building turns out inevitable.
Contrary to the previous studies, the approach developed here does not require instrumentation of the
building. Also, we consider buildings in normal occupation, that is in the presence of heating regulator,
which is an additional unknown of the problem. So, we identify a global system of the building together
with its heating system from:
1. meteorological data of the closest Meteo France station;
2. reference indoor temperature reconstructed by sector-based knowledge;
3. heating consumption collected either from Building Management Systems or from smart meters;
4. other heat supply (like lighting or presence of persons) estimated by sector-based and thermal
knowledge.
The identification is made at first by estimation of the parameters (7) defining the global model, by
minimizing one step prediction error. Then we adopted an approach of bayesian inversion which
provides a simulation of the posterior distributions of parameters and building indoor temperature.
The analysis of the stochastic simulations aims at studying the contribution of additional knowledge of
the problem (typical values of the parameters) and at showing the limits of the modeling hypotheses
in the case of certain real data sets.
Keywords: Identification, bayesian inversion, Gibbs sampler, stochastic simulation