Je Tiens

Keywords: Chauffe-eau Joule, Stockage d'energie, Electric hot water tank, Energy storage, Supply of hot water, Domestic water consumption, Dynamic Optimization, Multi-period model, MILP, Heuristics, Fokker-Planck

This thesis focuses on the development of advanced strategies for load shifting of large groups of electric hot water tanks (EHWT).

The first part of this thesis is dedicated to representing an EHWT as an input-output system. The idea is to design a simple, tractable and relatively accurate model that can be implemented inside a computing unit embedded in a "smart EHWT", for practical applications of optimization strategies. It includes in particular a phenomenological multiperiod model of the temperature profile in the tank and a model for domestic hot water consumption.

The second part focuses on the design of control strategies for a group of tanks. Three use-cases are studied. The first one deals with a small number of smart and controllable EHWT for which we propose a discrete-time optimal resolution method. The second use-case adresses a medium-scale group of controllable tanks for which we propose a heuristic to optimally schedule the heating periods. Finally, we present the modelling of the behavior of an infinite population of tanks under the form of a Fokker-Planck equation.

Résumé

Cette thèse s'intéresse au développement de stratégies de décalage de charge pouvant être appliquées à un parc de chauffe-eau Joule (CEJ).

On propose une modélisation entrée-sortie du système que constitue le CEJ. L'idée est de concevoir un modèle précis et peu coûteux numériquement, qui pourrait être intégré dans un "CEJ intelligent". On présente notamment un modèle phénoménologique multi-période d'évolution du profil de température dans le CEJ ainsi qu'un modèle de la demande en eau chaude.

On étudie des stratégies d'optimisation pour un parc de CEJ dont la résistance peut être pilotée par un gestionnaire central. Trois cas de figures sont étudiés. Le premier concerne un petit nombre de ballons intelligents et présente une méthode de résolution d'un problème d'optimisation en temps discret. Puis, on s'intéresse à un parc de taille moyenne. Une heuristique gardant indivisibles les périodes de chauffe (pour minimiser les aléas thermo-hydrauliques) est présentée. Enfin, un modèle de comportement d'un nombre infini de ballon est présenté sous la forme d'une équation de Fokker-Planck. As detailed in numerous studies, the increasing share of intermittent renewable electricity sources in the energy mix steadily complicates the management of electricity productionconsumption balance [START_REF]Energy roadmap[END_REF][START_REF] Edenhofer | IPCC special report on renewable energy sources and climate change mitigation[END_REF]. When used in addition to non-flexible means of production1 , these sources can even overload the grid, as witnessed during recent negative electricity price periods on French and German day-ahead markets [START_REF]Annual report[END_REF] 2 . Similar problems are observed at local levels in tension regulation across distribution grids originally designed for the sole purpose of electricity delivery under relatively steady conditions, which now have to put up photovoltaic production, causing further problems for electricity distribution companies.

Notations Notations relative to the EHWT

If the production of electricity is seldom flexible, then one may try to find some flexibility in the demand. This is the purpose of Demand Side Management (DSM), which is a collection of techniques aiming at modifying consumers' demand. DSM has an appealing potential [START_REF] Palensky | Demand side management: Demand response, intelligent energy systems, and smart loads[END_REF]. It appears as particularly relevant nowadays, since the global overcapacity of electricity production in Europe renders the construction of new flexible means of production non profitable, and therefore very unlikely.

A key factor in the growth of DSM is the availability of energy storage capacities. For this reason, network operators and electricity producers are searching and promoting new ways of storing energy. In this context, the large groups of electric hot water tanks (EHWT) found in homes in numerous countries 3 appear as very relevant, especially for the dominant problem of load-shifting applications. The main reasons supporting this fact are the intrinsic qualities of large storage: the large capacity of the population of EHWT, its geographically scattered characteristic, and its functioning.

In principle, heating of an EHWT can be freely scheduled, e.g. according to the price of electricity. While numerous advanced pricing policies have been studied and developed, the time-of-use pricing policy remains of dominant importance for most consumers, e.g. in France with the "night time switch"(which can be referred to as the historical strategy) 4 . Such policies define broad blocks of hours (for instance on-peak from 6:00 am to 10:00 pm, and off-peak from 10:00 pm to 6:00 am) during which a predetermined fixed rate is applied. The starting times of these periods being known, straightforward heating policies are commonly applied to each individual EHWT. Heating is turned-on immediately after reception of a wired communication signal broadcasted after the start of the off-peak period. Heating is turned-off when the EHWT is fully heated. This simple strategy increases electricity consumption in the night-time (one period when market electricity prices are low), while hot water is used in the next daytime. At large scales, the result is for the most part positive, but a negative effect is that the overall consumption of the group of EHWT rapidly decreases to a low level in the middle of the night, when the electricity production costs are the lowest, unfortunately.

This negative effect has to be addressed. With the fast-paced development of home automation, advanced heating strategies applied on large groups of EHWT are believed to enable further cost-reductions for both users and utilities [START_REF] Lane | A model of the domestic hot water load[END_REF]. The thesis presents works developed in this perspective, and aims at developing methods that can be key ingredients in the so-called "smart EHWT", which could integrate such advanced strategies.

Functioning and models of EHWT

A typical EHWT is a vertical cylindrical tank filled with water. A heating element 5 is plunged at the bottom end of the tank (see Fig. 1.1). The heating element is pole-shaped, and relatively lengthy, up to one third of the tank. Cold water is injected at the bottom while hot water is drained from the top at exactly the same flow-rate (under the assumption of pressure equilibrium in the water distribution system). Therefore, the tank is always full.

In the literature [Bla10, [START_REF] Kleinbach | Performance study of one-dimensional models for stratified thermal storage tanks[END_REF][START_REF] Zurigat | Influence of inlet geometry on mixing in thermocline thermal energy storage[END_REF], hot water storages are modeled as vertical columns driven by thermo-hydraulic phenomena: heat diffusion, buoyancy effects and induced convection and mixing, forced convection induced by draining and associated mixing, and heat losses at the walls. In the tank, layers of water with various temperature coexist (see Fig. 5.1 and Fig. 1.3). At rest, these layers are mixed only by heat diffusion which effects are relatively slow compared to the other phenomena [START_REF] Han | Thermal stratification within the water tank[END_REF].

The fact that of a non uniform temperature profile in the tank (increasing with height) remains in a quasi-equilibrium is called stratification [START_REF] Dincer | Thermal Energy Storage: Systems and Applications[END_REF][START_REF] Han | Thermal stratification within the water tank[END_REF][START_REF] Lavan | Experimental study of thermal stratified hot water storage tanks[END_REF]. In practice, this effect is beneficial for the user because hot water available for consumption is naturally stored near the outlet of the EHWT, while the rest of the tank is heated (see Fig. 1.2). The layer of high temperature gradient between cold and hot water (see Fig. 1.2) is called thermocline [START_REF] Zurigat | Influence of inlet geometry on mixing in thermocline thermal energy storage[END_REF]. The hot water consumption takes the form of a succession of drains of various magnitudes corresponding to the usages of the consumer (shower, bath, sanitation, cooking). The duration of the drains is short at the scale of the day, and their time of occurrence are not fixed but related to the inherent stochastic behavior of the user.

Most controlling strategies for groups of electric water storage do not account for the stratification and simply assume that the temperature in the tank is uniform [DL11, SCV + 13]. This approach is relevant for small-sized tanks, but are unable to make good use of the stratification phenomena. As will be shown in this thesis, including stratification modeling into those strategies can lead to non negligible performance improvements.

Early models of tanks including stratification have focused on large water storage tanks used in building basements for heating or air conditioning purposes [START_REF] Oppel | Computer simulation of stratified heat storage[END_REF]. Those tanks do not include heating (respectively chilling) elements, but are simply used to store large quantities of water heated (respectively chilled) by other means. Most authors use the cylindrical symmetry of the system and the fact that water flow is mostly vertical (due to the geometry of the system). The water temperature in the tank is assumed to be homogeneous at each height of the tank and the study is limited to one-dimensional models of various nature: convection-diffusion partial differential equations [START_REF] Oppel | Computer simulation of stratified heat storage[END_REF][START_REF] Zurigat | Stratified thermal storage tank inlet mixing characterization[END_REF], layer models [START_REF] Han | Thermal stratification within the water tank[END_REF], plug-flow models [START_REF] Kleinbach | Performance study of one-dimensional models for stratified thermal storage tanks[END_REF]. More recent work encompass heating components, such as heat exchangers in solar or thermodynamical water tanks [START_REF] Spur | Performances of modern domestic hot-water stores[END_REF][START_REF] Blandin | Modélisation et validation expérimentale de nouveaux concepts de ballons solaires à forte stratification[END_REF]. Three or two-dimensional (using rotational symmetry) models, often discretized for numerical simulation purposes such as computational fluid dynamics [Bla10, [START_REF] Han | Thermal stratification within the water tank[END_REF][START_REF] Johannes | Comparison of solar water tank storage modelling solutions[END_REF] or so-called zonal models based on the software TRNSYS [JFAR05, KBM10] can be found. These models, although accurate, are numerically intensive and mostly used in optimal design of exchangers and pipes found in most recent tanks.

A complexity trade-off must be found to reproduce the physical phenomena whose effects are observed in practice, while enabling fast determination of optimized strategies. With this aim in view, one-dimensional models appear as the ideal level of complexity. The works of modeling and control in this thesis will be limited to this case.

Control problems for EHWT

In principle, an EHWT can be seen as a two inputs, single output dynamical system (see Fig. 1.4). The two inputs are i) the heating power (which can be seen as a control variable) and ii) water outflow (or drain) chosen by the user. The internal state is the distribution of temperature of water in the tank, which can be used to define, for optimization purposes, performance indexes (outputs). The output can be loosely defined as the state of charge of the tank. As will appear in this thesis, a more detailed definition is required.

From an optimization viewpoint, the satisfaction of the user will be seen as a constraint, Heating power u(t)

Water drain which prevails over electric load management. In our case, the constraints can be defined using a comfort temperature set by the user. Water over this temperature can be blended with cold water, while water under the comfort temperature is useless. Additionally, functioning constraints can be considered. For instance, to prevent skin burns and EHWT malfunctions, a maximum (safety) temperature can be defined.

State of charge EHWT

The optimization problem can be defined from the electricity producer's perspective, by assuming that we have control of the power injection u(t) = (u 1 (t), ..., u k (t)) of a group of k tanks during a time period [t 0 , t f ]. The electricity producer usually desires to minimize a given objective function, while ensuring user's needs. For a given tank i, the fact that a strategy u i ensures user's comfort will be noted u i ∈ U i , U i being the set of admissible controls for tank i. A most useful problem of minimization concerns the cost of heating. Given a price signal for electricity over time c(t, •), an optimal control problem can be defined as follows:

min (u 1 ,...,u k )∈U 1 ×...×U k t f t=t 0 c(t, k j=1 u j (t)) k j=1 u j (t).
(P 1 )

We call this problem P 1 . In this formulation, the price depends on the total demand, since the whole population of EHWT is not negligible in the whole demand 6 . Alternatively, the objective function can be defined as a quadratic distance to an objective for the aggregated consumption f o (t). This gives min

(u 1 ,...,u k )∈U 1 ×...×U k t f t=t 0 (f o (t) - k j=1 u j (t)) 2 . (P 2 )
We call this problem P 2 . Works from the literature on these problems (or closely-related ones) generally consider that the temperature in the tank is uniform and neglect the stratification ([DL11, SCV + 13]). This assumption can lead to a violation of the comfort constraints caused by the supply of cold or tepid water, and do not take advantage of the stratification. For this reason, we wish to have a new look at these problems, and propose new solutions for them.

6 if the considered group is small, prices appear as flat and this problem splits into k distinct problems (since the functioning constraints are readily decoupled) min

u i ∈U i t f t=t0 c(t)u i (t).

Contributions of the thesis

The thesis proposes several contributions.

A first part is dedicated to representing an EHWT as an input-output system (see Fig. 1.4). The idea is to develop a simple, numerically tractable and relatively accurate model that can be implemented inside a computing unit embedded in a "smart EHWT" for practical applications of optimization strategies. In Chapter 2, a physics-based partial differential model including a natural convection term for the temperature is developed. Its terms are defined based on our physical understanding of the 1D dynamics in the tank. Fed with experimental recordings of data of electric consumption and water drains, the model is able to reconstruct the temperature profile in the tank. Then, this model is simplified in Chapter 3, into a phenomenological multi-period model which distinguishes several functioning states. This second model is computationally cheap and accurate. It may be integrated into the embedded chip of a smart EHWT. In the perspective of control design, a realistic domestic hot water consumption model is developed in Chapter 4. It can be used to randomly generate drain scenarios in near future, correlated with recent drains event. These scenarios can be used to proof-test optimal heating strategies which should avoid the prejudicial case of shortage of hot water. The temperature profile gives a precise description of the state of the EHWT, but is too heavy to be manipulated for control purposes. For this reason, we define three variables of interest that ease the handling of comfort constraints. This is the purpose of Chapter 5. This final modeling stage completes the definition of an input-output model for the EHWT subjected to a random hot water demand.

The second part focuses on the design of optimal control strategies for a group of tanks. It is assumed that the individual heating elements of the tanks can each be remotely piloted by a decision center (they are said to be controllable). Among them, upon request, a fraction can compute and transmit to the decision center information on their state. They are said to be smart. Three use-cases are studied. Chapter 6 adresses the case of a small number of smart and controllable EHWT (from 1 to 4). A discrete-time optimal resolution method for problems P 1 and P 2 is defined. In the case of a single unit, the optimization algorithm can be integrated in the chip of the tank, which will automatically define its heating strategy against an external price signal or with self-consumption objectives. In Chapter 7, we focus on a medium-scale group of controllable tanks (from a few hundreds to several millions) and present a heuristic for problem P 2 which keeps the heating period undivided to minimize thermo-hydraulic hazards. This heuristic, based on the stochastic nature of the problem, uses the statistical smoothing induced by the large number of tanks to approach the objective curve. If a fraction of the group is constituted of smart tanks, results are improved and can reach almost the same performance as in the case when all the tanks are smart (in which sub-optimality is less than 1%). Finally, the modeling of the behavior of a infinite population of tanks is presented in Chapter 8. This (prospective) approach aims at studying the output response of the whole group subjected to global strategies.

Perspectives and conclusions are discussed in Chapter 9.

Note. The works presented in this thesis have been the subject of the following publications and patents.

Introduction La gestion de la demande

Comme décrit dans de nombreuses études, la part grandissante des énergies renouvelables dans le mix énergétique complique la gestion de l'équilibre offre-demande de l'électricité [START_REF]Energy roadmap[END_REF][START_REF] Edenhofer | IPCC special report on renewable energy sources and climate change mitigation[END_REF]. Les sources renouvelabes peuvent s'ajouter à des moyens de production non-flexibles7 , et provoquer des surcharges du réseau comme en témoignent les épisodes récents de prix négatifs sur les marchés J+1 en France et en Allemagne [START_REF]Annual report[END_REF]. Des problématiques similaires sont observées au niveau local dans la régulation de tension du réseau de distribution, initialement dimensionné pour acheminer de manière relativement stable l'électricité au consommateur, qui doit désormais supporter une production photovoltaïque distribuée.

Pour compenser ce déficit de flexibilité sur les moyens de production d'électricité, on peut chercher d'autres gisements de flexibilité, notamment du cÃťté de la demande. C'est l'objet de la gestion de la demande (GD), un ensemble de techniques visant a modifier la demande des consommateurs. Le potentiel de la GD apparaît prometteur [START_REF] Palensky | Demand side management: Demand response, intelligent energy systems, and smart loads[END_REF]. C'est d'autant plus vrai actuellement, car la situation de surcapacité de la production d'électricité en Europe rend la construction de moyens de production flexibles non rentable et par conséquent peu probable dans les années à venir.

Un facteur clé pour le développement de la GD est la disponibilité de capacités de stockage. Pour cette raison, les gestionnaire de réseaux et les producteurs d'électricité recherchent de nouveaux moyens de stockage. Dans ce contexte, les importants parcs de chauffe-eau Joule (CEJ) domestiques de nombreux pays8 apparaissent prometteurs, particulièrement pour le problème dominant du décalage de charge. Les principales raisons supportant ce fait sont les qualités intrinsèques de ce stockage : la grande capacité des parcs de CEJ, leur caractère réparti sur le territoire, et leur fonctionnement.

En principe, la chauffe d'un CEJ peut être planifiée librement, par exemple en fonction des prix de l'électricité. Si de nombreuses politiques tarifaires ont été étudiées et élaborés, la tarification en fonction de l'heure de consommation reste dominante pour la plupart des consommateurs, par exemple en France avec le système d'heures pleines/heures creuses 9 . Ces politiques définissent des plages horaires (par exemple des heures pleines de 6h à 22h, et des heures creuses de 22h à 6h) pendant lesquelles un tarif prédéterminé est appliqué. Le début de ces plages étant connu, des stratégies de chauffe simples sont généralement appliquées à chaque CEJ. La chauffe est mise en marche immédiatement après réception d'un signal analogique correspondant à la plage d'heures creuses. La chauffe est ensuite arrêtée lorsque toute l'eau du CEJ est chaude. Cette stratégie élémentaire favorise naturellement la consommation d'électricité pendant la nuit (une période pendant laquelle les prix de marché de l'électricité sont bas), alors que l'eau chauffe est utilisée la journée suivante. A grande échelle, le résultat est dans l'ensemble positif, mais un effet indésirable est que la consommation globale redescend rapidement à un bas niveau au milieu de la nuit, lorsque les coûts de production de l'électricité sont les plus bas.

Cet effet négatif doit être corrigé. Le développement rapide de la domotique ouvre la voie à des stratégies de chauffe nouvelles sur de grands parcs de CEJ, qui peuvent s'avérer avantageuses à la fois pour les consommateurs et pour les producteurs d'électricité [START_REF] Lane | A model of the domestic hot water load[END_REF]. Cette thèse présente des travaux développés dans cette perspective, et vise à développer des méthodes qui pourraient être intégrées à des "CEJ intelligents" mettant en oeuvre de telles stratégies.

Fonctionnement et modèles de CEJ

Un CEJ prend usuellement la forme d'un ballon cylindrique vertical, rempli d'eau. Une résistance chauffante 10 est plongée dans le bas du ballon (voir Fig. 1.5).

La résistance est longiligne et relativement et relativement grande (elle couvre jusqu'à un tiers de la longueur du ballon). Lorsque l'utilisateur soutire de l'eau, cette eau est prélevée dans le haut du ballon pendant que de l'eau froide est injectée dans le bas du ballon au même débit (sous l'effet de l'équilibre de pression du réseau de distribution d'eau), de telle sorte que le ballon est toujours plein.

Dans la littérature [Bla10, KBK93, ZLG91], les réservoirs d'eau chaude sont modélisés comme des colonnes verticales, régi par des phénomènes thermo-hydrauliques : diffusion thermique, effets de flottabilité dus à la variation de densité de l'eau en fonction de sa température et des mélanges qu'ils engendrent, convection due au soutirage et mélange engendré, et enfin pertes thermiques vers l'extérieur. Dans le ballon, des couches d'eau de températures différentes coexistent (voir Fig. 5.1 et Fig. 1.7). Au repos, ces couches ne sont affectées que par la diffusion thermique, dont les effets sont relativement lents comparés aux autres phénomènes [START_REF] Han | Thermal stratification within the water tank[END_REF].

Le fait qu'un profil (croissant) de température dans le ballon, non uniforme, reste dans un quasi-équilibre est appelé phénomène de stratification [DR10, HWD09, LT77]. En pratique, cet effet est utile pour l'utilisateur, car l'eau chaude disponible pour la consommation est naturellement stockée près de l'évacuation du CEJ, pendant que le reste du ballon est chauffé (voir Fig. 1.6). La zone de fort gradient de température entre l'eau chaude et l'eau froide (voir Fig. 1.6) est appelée thermocline [START_REF] Zurigat | Influence of inlet geometry on mixing in thermocline thermal energy storage[END_REF]. En pratique, la consommation d'eau chaude prend la forme d'une succession de soutirages d'amplitudes variables, correspondant aux usages du consommateur (douche, bain, hygiène, cuisine). La durée de ces soutirages est courte à l'échelle de la journée, et les moments où ils ont lieu ne sont pas fixés, mais inhérents au caractère stochastique du comportement de l'utilisateur.

La plupart des stratégies de contrôle de parcs de CEJ ne prennent pas en compte la stratification et supposent simplement que la température dans le ballon est uniforme [DL11, SCV + 13]. Cette approche est pertinente pour les ballons de petites tailles, mais ne permet pas de valoriser le phénomène de stratification. Comme il sera exposé dans cette thèse, inclure un modèle de stratification dans ces stratégies peut permettre des améliorations Dans la litt Ãľrature, les premiers modèles de ballon incluant la stratification se sont intéress Ãľs à des réservoirs de grande taille dans les sous-sols des bâtiments utilisés pour le chauffage ou la climatisation [START_REF] Oppel | Computer simulation of stratified heat storage[END_REF]. Ces réservoirs n'incluent pas de résistance chauffante ou de système de réfrigération direct. Ils sont simplement utilisés pour stocker de grandes quantités d'eau chauffée ou refroidie par d'autres moyens. La plupart des auteurs utilisent la symétrie de révolution du système et le fait que l'écoulement est principalement vertical (en raison de la géométrie du système). La température du ballon est alors supposée homogène pour une hauteur donnée, et l'étude est limitée à des modèles unidimensionnels de différentes natures : équations aux dérivées partielles de convectiondiffusion [START_REF] Oppel | Computer simulation of stratified heat storage[END_REF][START_REF] Zurigat | Stratified thermal storage tank inlet mixing characterization[END_REF], modèles par couches [START_REF] Han | Thermal stratification within the water tank[END_REF], modèles piston [START_REF] Kleinbach | Performance study of one-dimensional models for stratified thermal storage tanks[END_REF]. Les travaux les plus récents comprennent des élements chauffants, comme des échangeurs thermiques dans les chauffe-eau solaires ou thermodynamiques [START_REF] Spur | Performances of modern domestic hot-water stores[END_REF][START_REF] Blandin | Modélisation et validation expérimentale de nouveaux concepts de ballons solaires à forte stratification[END_REF]. Des modèles tri-ou bidimensionnels (en utilisant la symétrie de révolution), souvent discrétisés pour des objectifs de simulations numérique (comme par exemple de mécanique des fluides numérique) [Bla10, HWD09, JFAR05], ou des modèles dits zonaux basés sur le logiciel TRNSYS [JFAR05, KBM10] sont l'objet de certaines Ãľtudes. Ces modèles, bien que précis, sont coûteux numériquement et principalement utilisés dans la conception optimale (hors ligne) des échangeurs et des tuyaux dans les chauffe-eau les plus récents.

Un arbitrage entre précision et complexité doit être établi pour reproduire les phénomènes physiques observés en pratique, tout en permettant le calcul rapide de stratégies optimisées. Dans cette optique, les modèles unidimensionnels apparaissent comme le niveau idéal de complexité. Les travaux de modélisation et de contrôle de cette thèse se limiteront à ce cas.

Problèmes de contrôle pour les parcs de CEJ

En principe, un CEJ peut être vu comme un système dynamique a deux entrées et une sortie (voir Fig. 1.8). Les deux entrées sont i) la puissance de chauffe (qui est une variable de contrôle) et ii) le débit d'eau sortant soutiré par l'utilisateur. L'état interne est le profil de température de l'eau dans le ballon, qui peut être utilisé pour définir, dans un but d'optimisation, des indices de performances (sorties). Ces sorties peuvent (pour l'instant approximativement ) être définies comme l'état de charge du CEJ. Comme on le verra dans cette thèse, une définition plus précise est nécessaire.

Du point de vue de l'optimisation, la satisfaction de l'utilisateur est une contrainte, qui prévaut sur la gestion de la charge électrique. Dans notre cas, les contraintes peuvent être définies en utilisant une température de confort réglée par l'utilisateur. L'eau au dessus de cette température peut être mitigée, en la mélangeant avec de l'eau froide, alors que l'eau en dessous de la température de confort est inutile. De plus, des contraintes de fonctionnement sont a prendre en compte. Par exemple, pour éviter des brûlures cutanées ou des dysfonctionnements du CEJ, une température maximum (de sécurité) peut-être définie.

Le problème d'optimisation peut être défini du point de vue du producteur, en supposant qu'on dispose du contrôle de l'injection de puissance u(t) = (u 1 (t), ..., u k (t)) d'un parc de k ballons pendant une période [t 0 , t f ]. Le producteur d'électricité souhaite généralement minimiser une certaine fonction objectif, toute en satisfaisant les contraintes de l'utilisateur. Pour un ballon donné i, le fait qu'une stratégie u i garantisse le confort de l'utilisateur sera noté u i ∈ U i , U i étant l'ensemble des contrôles admissibles pour le ballon i. Un problème important de minimisation concerne le coût de la chauffe. Étant donné un signal prix pour l'électricité au cours du temps c(t, •), un problème de contrôle optimal peut être défini comme min

(u 1 ,...,u k )∈U 1 ×...×U k t f t=t 0 c(t, k j=1 u j (t)) k j=1 u j (t).
(P 1 )

On appelle ce problème P 1 . Dans cette formulation, le prix dépend de la consommation totale en électricité, dans la mesure ou la taille des parcs de CEJ n'est pas négligeable dans la demande globale 11 .

11 si le groupe considéré est de petite taille, étant donné que les contraintes de fonctionnement sont découplées, les prix apparaissent comme plats et ce problème se décompose en k problèmes distincts min

u i ∈U i t f t=t0 c(t)u i (t).
Alternativement, la fonction objectif peut-être définie comme une distance quadratique à un objectif pour la consommation agrégée f o (t). Le problème est alors min

(u 1 ,...,u k )∈U 1 ×...×U k t f t=t 0 (f o (t) - k j=1 u j (t

Contributions de la thèse

Cette thèse propose plusieurs contributions. Une première partie est dédiée à la représentation d'un CEJ en tant que système entrée-sortie (voir Fig. 1.8). L'idée est de développer un modèle simple, adaptable et assez précis, qui puisse être utilisé dans une unité de calcul embarquée dans un CEJ intelligent pour l'application pratique des stratégies d'optimisations. Dans le chapitre 2, un modèle physique prenant la forme d'un système d'équations aux dérivées partielles est développé, incluant notamment un terme de convection naturelle. Les termes du modèles sont basés sur notre compréhension physique de la dynamique 1D dans le ballon. Validé par des données expérimentales liées à des chroniques de soutirage et de consommation électrique, le modèle se montre capable de reconstruire le profil de température Ãă l'intérieur d'un ballon. Dans un second temps, ce modèle est simplifié dans le chapitre 3 en un modèle phénoménologique multi-périodes qui distingue plusieurs états de fonctionnement. Ce second modèle est précis et peu coûteux numériquement. Il est possible de l'intégrer dans un calculateur embarqué dans un CEJ intelligent. Dans la perspective de la conception des stratégies de contrôle, un modèle réaliste de consommation d'eau chaude sanitaire est développé dans le chapitre 4. Il peut être utilisé pour générer aléatoirement des scénarios de soutirages dans le futur proche, corrélé avec l'historique récent de soutirage. Ces scénarios peuvent être utilisés pour mettre à l'épreuve des stratégies de chauffe optimale et vérifier qu'elles ne génèrent pas une pénurie d'eau chaude. Le profil de température offre une description précise de l'état du CEJ, mais son caractère distribué le rend trop lourd pour être manipulé dans un but de contrôle. Pour cette raison, on définit trois variables d'intérêt qui facilitent la manipulation des contraintes de confort. C'est l'objet du chapitre 5. Cette étape finale de modélisation termine la définition d'un modèle entrée-sortie pour un CEJ soumis à une demande d'eau chaude sanitaire aléatoire.

La seconde partie de la thèse s'intéresse à la conception de stratégies de contrôle optimal pour un parc de CEJ. On suppose que la chauffe de chaque ballon peut être pilotée à distance par un centre de décision (les ballons étant dits contrôlables). Parmi eux, une partie peut calculer et transmettre sur demande au centre de décision des informations sur leurs états. Ces ballons sont dits intelligents. Trois cas de figure sont étudiés. Le chapitre 6 s'intéresse aux cas d'un petit nombre de ballons contrôlables et intelligents (de 1 à 4). Une méthode de résolution optimale en temps discret pour P 1 et P 2 y est définie. Dans le cas d'un ballon unique, l'algorithme d'optimisation peut être intégré dans le calculateur embarqué du ballon, et définit automatiquement sa stratégie de chauffe future en utilisant par exemple un signal prix extérieur ou des objectifs d'auto-consommation. Dans le chapitre 7, on étudie des parcs de tailles intermédiaires (de quelles centaines à plusieurs millions). On propose une heuristique de résolution du problème P 2 , gardant les périodes de chauffes indivisées pour minimiser les aléas thermo-hydrauliques. Cette heuristique repose sur la nature stochastique du problème, et utilise le lissage statistique naturel généré par le grand nombre de ballons pour se rapprocher de la courbe objectif. Si une partie des ballons est contrôlable et intelligente, les résultats s'améliorent et la perte d'optimalité peut descendre à moins de 1%. Enfin, un modèle de comportement d'un parc infini de ballons est proposé dans le chapitre 8. Cette approche (prospective) vise à étudier la réaction d'un parc soumis à des stratégies d'ensemble.

Le chapitre 9 est consacré à la conclusion, et à la présentation de perspectives.

Note. Les travaux présentés dans cette thèse ont été l'objet des publications et brevets suivants. In this chapter, we develop a distributed parameter model for the dynamics of the temperature profile in an EHWT. The microscopic and macroscopic effects observed in response to the water drain and the power injected in the tank take central roles in this model.

Liste des papiers de conférence

The model can be seen as an extension of existing one-dimensional convectiondiffusion linear equations modeling the draining convection and mixing originally developed in [START_REF] Oppel | Computer simulation of stratified heat storage[END_REF][START_REF] Zurigat | Stratified thermal storage tank inlet mixing characterization[END_REF][START_REF] Zurigat | Influence of inlet geometry on mixing in thermocline thermal energy storage[END_REF] in the general context of storage tanks. In details, to the classic governing equation already found in the previously cited works, we add a nonlinear velocity term stemming from an empirical law representing turbulent natural convection caused by heating, and we explicitly include heating power as a source term1 .

The model we develop is shown to be in accordance with experimental data presented in this chapter. These data clearly stress the following: i) heating water with an EHWT takes time, ii) during the heating process, a zone of uniform temperature appears and grows until it covers the whole tank, iii) draining induces a piston flow, and also causes some internal mixing which is non negligible. As will be shown, the proposed distributed parameter model is able to reproduce these observations. The chapter is organized as follows. After having described the proposed model in Section 2.1, we illustrate it by means of simulations and compare it against experimental data in Section 2.2. A summary is given in Section 2.3.

PDE for heating and draining

Draining model as a PDE

To emphasize the effects of stratification, our model solely uses one-dimensional partial differential equations. The works of Zurigat on draining effects in stratified thermal storage tanks [START_REF] Zurigat | Stratified thermal storage tank inlet mixing characterization[END_REF][START_REF] Zurigat | Influence of inlet geometry on mixing in thermocline thermal energy storage[END_REF] serves as baseline. The novelty is to introduce the heating system. It is treated in § 2.1.2 and § 2.1.3. The equation below accounts for draining and its induced turbulent mixing effects. It is an usual one-dimensional energy balance where the turbulence is lumped into a diffusion term

∂ t T + ∂ x (v d T ) = (α th + α d )∂ xx T .
In this equation, T (x, t) is the temperature at time t and height x, v d ≥ 0 is the velocity induced by the draining (assumed to be spatially uniform but time-varying), α th is the thermal diffusivity and α d is an additional turbulent diffusivity term representing the mixing effects. Zurigat [START_REF] Zurigat | Stratified thermal storage tank inlet mixing characterization[END_REF][START_REF] Zurigat | Influence of inlet geometry on mixing in thermocline thermal energy storage[END_REF] considers the same equation and introduces the ratio d = α th +α d α th . An experimental correlation is shown with the Reynolds number and the Richardson number Ri in the tank2 defined as

Ri = gα V (T in -T a )L c v 2 d (2.1)
where g is the gravitational acceleration, α V is the volumetric coefficient of thermal expansion of the fluid (here water), T in is the temperature of the inlet water, T a is the ambient temperature and L c is a characteristic vertical length3 . This correlation can be used in our case.

A heat losses term (to the exterior of the tank assumed to be at temperature T a ) can be added to this equation. Then, one obtains

∂ t T + ∂ x (v d T ) = d α th ∂ xx T -k(T -T a ) (2.2)
where the factor k is defined by

k = U h Π
Sρc p where ρ is the density of water, c p its specific heat capacity, U h is the overall heat transfer coefficient based on the tank internal surface area, Π is the tank internal perimeter and S its effective cross-section.

Note h the vertical length of the tank internal volume. Equation (2.2) is assumed to hold over Ω × I, where

I =]t 0 , t f ] is a time interval and Ω =]0, h[. Classically, we consider boundary conditions of the Robin [DL93] form v d (T (0, t) -T in ) + d α th ∂ x T (0, t) = 0 at x = 0, and of the Neumann [DL93] form ∂ x T (h, t) = 0 in x = h,
meaning that energy is allowed to leave the system with the outlet flow but not with diffusion.

Including heating and buoyancy forces

Equation (2.2) integrates the most obvious phenomena taking place in an EHWT. The effects of heating are of two types: direct and indirect. The direct effect is a source term in the balance equation. The indirect effect is buoyancy. It can be added in the model at the expense of linearity, as is described below.

When the heating system is on, temperature of water around the heating element starts to rise. By buoyancy, hot water replaces colder water above by a phenomenon called Rayleigh-Bénard convection [START_REF] Peterson | Thermal Aspects of Nuclear Reactors[END_REF][START_REF] Gauthier | Convection turbulente dans une cellule de Rayleigh-Bénard cryogénique : de nouveaux éléments en faveur du Régime Ultime de Kraichnan[END_REF][START_REF] Gibert | Convection thermique turbulente : Panaches et fluctuations[END_REF]. This convection can take various forms depending on the characteristics of the system, represented by the Rayleigh number

Ra = gα V (T s -T ∞ )L 3 c να th
where ν is the kinematic viscosity of the fluid, T s is the temperature of the surface (here the heating element), and T ∞ is the temperature in the tank far from it. This adimensional number scales the effects of buoyancy and conduction: if it is low, the conduction will be the main heat transfer factor, if it is high, the natural convection will predominate. Over a critical value (Ra = 1108 for [START_REF] Peterson | Thermal Aspects of Nuclear Reactors[END_REF], Ra = 3, 5 • 10 4 for [START_REF] Gauthier | Convection turbulente dans une cellule de Rayleigh-Bénard cryogénique : de nouveaux éléments en faveur du Régime Ultime de Kraichnan[END_REF]), turbulent natural convection appears under the form of a pattern of plumes forming convection cells called Bénard cells which can take various forms and sizes. For any EHWT found in households, even with a small T s -T ∞ difference, the Rayleigh number is far over the critical value, and plumes of turbulent water appears over the heating system (see Fig. 2.1 reproduced from [START_REF] Blandin | Modélisation et validation expérimentale de nouveaux concepts de ballons solaires à forte stratification[END_REF]). Therefore, convection dominates conduction. To include this effect into (2.2), we simply consider that, at each given height, two distinct temperatures co-exist in the convection cells. Then, to our equation on T , we append an interacting equation bearing on a new physical quantity ∆T (x, t) representing the temperature spread at each height x over T . This gives the following system

∂ t T + ∂ x (v d T ) = d α th ∂ xx T + Φ∆T -k(T -T a )
(2.3)

∂ t ∆T + ∂ x ((v d + v nc )∆T ) = d α th ∂ xx ∆T -Φ∆T + P W .
(2.4)

In the two equations above, three terms have been added: a velocity term v nc of natural convection, which is responsible for transport of energy in the system, a heat exchange term Φ(x, t)∆T (x, t) (representing at each height the mixing induced by natural convection being proportional to the temperature spread), and the spatially distributed source term P W (representing the power injected in the tank all along the element, generally from the bottom of the tank to one third of its height), which drives the dynamic of ∆T . The boundary conditions for (2.3) remain unchanged, while the boundary conditions for (2.4), as a temperature spread, are v d ∆T (0, t)

+ d α th ∂ x (∆T )(0, t) = 0 at x = 0 and ∂ x (∆T )(h, t) = 0 at x = h.

Model for natural convection and internal heat transfer

We now introduce a model for the transport velocity appearing in (2.4). This velocity v nc (x, t) is non-constant. It is non-zero at a given altitude x only if there exists colder water over the height x (i.e. downstream). We give to v nc the following integral form

v nc (x, t) = v( h x [T (x, t) + ∆T (y, t) -T (y, t)] ζ + dy) ξ (2.5)
where [z] + is the positive part of z and v is a positive factor, and where 0 < ζ < 1 and 0 < ξ < 1 are tuning parameters the value of which are chosen to fit experimental data. These parameters reduce the impact of the downstream temperatures differences and smooth the velocity when it is nonzero. The exchange coefficient Φ between the two equations is also non-constant and we model it as

Φ(x, t) = φ[v max -v nc (x, t)] + (2.6)
where φ and v max are two tuning parameters. The rationale behind this expression is that the horizontal mixing is stronger when the natural convection flow reaches the upper part of the tank and has a lower speed.

Summary of the model

According to the previous discussion, the EHWT can be represented by two distributed state variables, T and ∆T , governed by (2.3) and (2.4). In those governing equations two velocities appear: v d which is spatially uniform and is equal to the output flowrate (drain) of the system, and v nc which is defined in (2.5) to model the effects of natural convection. Heat is injected into the system through a distributed source term P W and the heat exchange between the two equations is proportional to ∆T with a coefficient (variable in space) defined in (2.6). Finally, α th , d , k, v, ζ, ξ, φ and v max are constant parameters depending on physical constants and the geometry of the tank. Typical values for the EHWT defined in Table 2.1 are given in Table 2.2. They result from an identification procedure.

Model validation

Experimental setup

α th m 2 s -1 1.43 • 10 -7 d - 13 k s -1 1.43 • 10 -6 v m 1-ξ •s -1 K -ζξ 10 -3 ζ - 0.2 ξ - 0.5 φ m -1 0.03 v max ms -1 0.35 Table 2
.2: Parameters of the model for the associated EHWT.

The water tank has been equipped with internal temperature sensors recording temperature at 15 locations of different heights, 15 cm deep inside the water tank (see Fig. 1.1). This depth is sufficient to bypass the insulation of the tank. It is assumed that the sensors have no effect on the flows (e.g. that they do not induce significant drag). Besides, the following quantities have been recorded with external sensors: injected power, water flow at the inlet, water temperature at the inlet. These three quantities feed the model, the output of which can be compared with the temperature measured by the sensors. The comparisons are directed into an optimization procedure identifying the coefficients given in Table 2.2. Conducted experiments took the form of fourteen 24 h runs with a sampling rate of 1Hz. Histories for drain are taken from the normative sheets emitted by the French norm organism [NF 11] for a tank of such capacity, associated with a classical night-time heating policy until total load. Subsequent experiments consider similar total consumption but with different drain/heat combinations and overlaps to test the model under various situations.

Results

For sake of illustration, several operating conditions are reported next. Simulations have been conducted on a quad-core Intel Core i7-4712HQ processor equipped with 16Go of memory. Numerically, this system of equations can be solved with finite difference schemes. We use a Crank-Nicholson scheme (see [START_REF] Allaire | Numerical Analysis and Optimization[END_REF]) for the diffusive term, and an upwind scheme for the convective part. The later is stable only conditionally to a Courant-Friedrichs-Lewy condition [START_REF] Allaire | Numerical Analysis and Optimization[END_REF]. In our case, due to the non-linear nature of v nc , short time-steps must be chosen. In turn, this increases the computational load which is already high due to the evaluation of the integral appearing in v nc , for each space-step. This can lead to long computation times (see Table 2.3). Fig. 2.2 (a) shows the variation of the temperature during a heating period for a tank initially completely cold. Fig. 2.2 (b) shows the response to draining of a heated tank. The results are quite satisfactory. Some mismatch appears in the lower part of the tank during draining, and in upper parts during heating. They can be due to the chosen one-dimensional representation, since the neglected effects of radial inhomogeneity may be stronger near the ends of the tank, and to some mixing effects that have not be taken into account. However, the results show that the model is accurate enough, even in 24 h open-loop runs. To support this statement, the distribution of the absolute difference between experimental value and model prediction is given in Table 2.3 (produced over the whole set of data).

Err. 0-2 o C 2-4 o C 4-6 o C 6-8 o C

Summary

In this chapter, we have presented a first model for an EHWT with consumed water and injected power as inputs. This model accounts for the direct and indirect effect of heating. Comparisons with experimental data reveal that the model is capable of reproducing with a relative degree of accuracy the transient behaviors. A closer look at both numerical simulations and experimental results highlight the formation of a spatially uniform temperature distribution which gradually extends itself upwards to the top of the tank (see Fig. 2.2 (b)). We consider that buoyancy induced forces, by generating a local natural convection phenomenon, are the root cause for this. The homogeneous zone is followed by an increasing profile of temperature in the upper part of the tank, staying untouched due to stratification if heat diffusion is neglected. This property suggests an alternative model focused on the dynamics of the uniform temperature of the zone with an ordinary differential equation (ODE). This simplification is presented next, in Chapter 3.

Chapter 3

Multi-period dynamical modeling for electric hot water tank

Modèlisation dynamique multi-periode des CEJ. Dans ce chapitre, nous simplifions le modèle précédent en distinguant trois régimes de fonctionnement. Pour chaque régime, la dynamique est décrite sous la forme d'une equation aux dérivées partielles ou aux dérivées ordinaires. La présentation du modèle est suivie d'une validation expérimentale.

The physics-based model presented earlier in Chapter 2 is concise and relatively satisfactory. However, the accuracy, and, above all, the computational cost associated to its numerical resolution can be seen as strong limitations if one desires to embed it in a "smart EHWT". In this chapter, we decouple heating and draining effects and develop a new model, based on the decomposition of the dynamics according to the dominant effect at stake. We distinguish three phases (periods): heating, draining and rest.

For heating, we reproduce the behavior observed earlier in the experimental data in which the temperature increases first at the bottom of the tank forming a spatially uniform temperature distribution which gradually extends itself upwards to the top of the tank. This homogeneous zone is followed by an increasing profile of temperature in the upper part of the tank, remaining untouched due to stratification (heat diffusion being neglected in this case). As previously, draining is treated as a convection parameter and its associated mixing effects are reproduced by a diffusion term. However, combined to the convection-diffusion equation, we model the effects of the water nozzle which creates a mixing zone of varying temperature and volume. The cascade represents a typical Stefan problem [START_REF] Fasano | General free-boundary problems for the heat equation[END_REF][START_REF] Fasano | General free-boundary problems for the heat equation, ii[END_REF]. Finally, rest phases are simply driven by diffusion and losses. Sequencing the three phases constitutes a multi-period model. This multi-period model is the main contribution of this chapter.

The chapter is organized as follows. Preliminary observations are displayed in Section 3.1. Section 3.2 is dedicated to the presentation of the multi-period model which is the main contribution of the chapter. Comparative studies reported in Section 3.3 conclude that this multi-period model is more accurate and more computationally efficient than the physics-based model presented in Chapter 2. A summary is given in Section 3.4.

Preliminary observations

The comparison of the physics-based model of § 2.1.4 against experimental data is overall satisfactory even in 24 h open-loop runs, but close-up inspections have revealed some possibilities of improvement. We now detail these.

Firstly, numerical results of this model and experimental data concur and clearly show, during heating, the creation of a temperature "plateau" starting from the bottom of the tank (see Fig. 3.4 (a.1)). In numerical simulations, this plateau has increasing temperature and length as heating goes on, but leaves temperatures at greater heights untouched and only progressively covers the whole of the tank. This phenomena is seen in experimental data in the exact same way at the exception of a small temperature backward flow observed in the highest region of plateau (see Fig. 3.4 (a.1)). These observations support the validity of the model, but they also suggest some simplifications could be performed.

Secondly, a mismatch appears during draining in the lower part of the tank. This mismatch consists of an underestimation of the injected water temperature and a shift of the location of the thermocline (see Fig. 3.4 (a.2)). It is believed that the mismatch arises from a local mixing in the bottom of the tank, such as the one coming from a strong squirt out of the water injection nozzle. In turn, the dynamics of such mixing and its effect on the thermocline location have a strong effect on the temperature profile that cannot be easily reproduced by the simple convection-diffusion equation in a fixed domain. Accounting for them calls for a transformation of the fixed domain into a time-dependent one and an adaptation of boundary conditions.

In practice, it appears that draining and heating effects on the tank are mostly taking place over disjoint periods. Therefore, the time interval during which the system is considered can be split into distinct subintervals, or periods. Introducing distinct dynamics for each period offers the double advantage of reducing the computational burden (enabling tailoring numerical schemes for each dynamics) and of bringing significantly improved flexibility compared to a single system of PDE.

We now present the multi-period model we propose.

Multi-period model for heating, draining and heat losses

Separations into non overlapping periods

Consider an initial time t 0 , an initial temperature profile T 0 (x), and a time t f at which one wants to determine the temperature profile T (•, t f ). On [t 0 , t f ], the tank is submitted to draining and heating, characterized by the draining velocity v d (t) and the injected power u(t) (related to the previously defined P W via the relation u(t) = h 0 P W (y, t)dy). We assume that v d and u are piecewise constant and left-continuous (at each discontinuity point).

Let us define T u = (t u 0 = t 0 , ..., t u mu = t f ) and T v = (t v 0 = t 0 , ..., t v mv = t f ) the sequence of discontinuity points respectively of u and v d , and T = T u ∪T v = (t 0 , ..., t m ) the sequence of discontinuity points of u and v d , such that t 0 < t 1 < ... < t m = t f . This sequence T defines a succession of m time intervals ]t i , t i+1 ] of length ∆t i . In each time interval, the tank is in one (and only one) of the phases I, II, IIIa, IIIb defined below. Over each time interval, u(t) and v d (t) are constant (as illustrated in Fig. 3.1). Over each interval ]t i , t i+1 ], we desire to determine the temperature profile as a function of time, in particular at final time t i+1 . The profile at the end of a phase serves as initial condition for the following phase. Three functions F I , F II , F III (accounting for both IIIa and IIIb) map an initial profile and working conditions to temperature profiles for future times. We note,

Draining

v d (t) = 0 v d (t) > 0 Heating u(t) = 0 I IIIa u(t) > 0 II IIIb Draining velocity v d (m • s -1 ) Heating power u (W ) Configuration Time Time Time t v 0 t 0 u max I I II II IIIa IIIa IIIb IIIb II t 1 t 2 t 3 t 4 t 5 t 6 t 7 t 8 t 9 = t f t v 1 t v 2 t v 3 t v 4 t v 5 t v 6 t v 8 = t f t u 1 t u 2 = t f t u 0 t v 7
• T (•, t i+1 ) = F I (T (•, t i ), ∆t i ) • T (•, t i+1 ) = F II (T (•, t i ), ∆t i , u) • T (•, t i+1 ) = F III (T (•, t i ), ∆t i , v d , u).
Clearly, if one wishes to compute the temperature profile at any time of interest t f , one only needs to compute the sequence of intermediate profiles T (•, t i ), i = 1, ..., m -1 as a function of the previous ones by a chain rule. Note that for the computation on a short interval [t 0 , t f ] (e.g. if one considers a succession of nearby times of interest), T can be reduced to a short list of events. Interestingly, a comparable split is developed in [START_REF] Kreuzinger | Mathematical modelling of a domestic heating system with stratified storage tank[END_REF] for the case of a water storage tank with external heating.

We now detail these mappings, for any index i.

Phase I: Rest

In this part, we consider periods without any draining or heating.

Physical considerations

The only phenomena driving the temperature profile are diffusion and heat losses.

Dynamics

The input variables of F I are the initial profile, say T 0 (•), and the duration ∆t i . They serve in the following diffusion-heat losses one-dimensional system:

∂ t T = α th ∂ xx T -k(T -T a ) on Ω × I ∂ x T (0, t) = 0 on I ∂ x T (h, t) = 0 on I T (x, 0) = T 0 (x) on Ω              (3.1)
where I =]0, ∆t] and Ω = [0, h].

We have

F I (T (•, t i ), ∆t i ) = T (•, t i+1 ).

Numerical considerations

Numerically, this system can be solved relatively easily with finite difference schemes. We use a Crank-Nicholson scheme [START_REF] Allaire | Numerical Analysis and Optimization[END_REF] on a linearly spaced mesh.

Phase II: Heating

Physical considerations

Heating modeling can be simplified thanks to the plateau discussed earlier. Turbulence generated by buoyancy effects during the heating process is the cause of a local mixing.

Here, we consider that this mixing is perfect on the plateau which is an area [0, x p (t)], and that the buoyancy effects do not affect stratification in heights above x p (t). For sake of simplicity of presentation, only the case without heat losses is exposed here, but losses can be included without too much difficulty (this is actually done for the simulation presented in Section 3.3). To simplify the dynamics, the diffusion phenomena have to be neglected. Then, the governing equations take the form of an ODE that we derive below.

Dynamics

The input variables of F II are the initial profile T 0 (•), the duration ∆t i , and the constant value u of the heating power. The plateau temperature is noted T p (t). It is related to x p (t) by the equation

T p (t) = T 0 (x p (t)) (3.2)
corresponding to the continuity assumption at the interface between the plateau and the (untouched) initial profile (see Fig. 3.2). An energy balance (illustrated in Fig. 3.2) gives

T p (t)x p (t) = xp(t) 0 T 0 (x)dx + t 0 u(s) Sρc p ds. ( 3
dx p dt = u Sρc p x p T 0 (x p ) , x p (0) = 0, t ∈]0, ∆t] (3.4)
which directly gives T p (t) using (3.2).

For completeness, other phenomena can be included:

• Heat losses at walls (with losses coefficient k) is a local phenomena that does not alter the shape of a plateau. Its effects are some decrease of T 0 (•) towards ambient temperature T a in the form of a exponential factor e -kt of the initial profile and a additional term in the dynamics of x p (t).

• If the plateau is not exactly uniform but a trend can be observed (for instance repeatable variations around the heating elements), one can separate the temperature of the plateau in two components with distinct arguments T p (x, t) = T pt (t) + T px (x) and then study the dynamics of T pt (t).

• Finally, the small backward energy flow that is always observed (see Fig. 3.4 (a.1)) at the interface (in a stronger way at the top of the tank) can be modeled by relaxing the continuity hypothesis (3.2) and replacing it with

T p (x p (t), t) + T p∆ (x p (t)) = T 0 (x p (t)) (3.5)
where the continuity gap T p∆ depends on the geometry of the tank (and has to be identified).

Integration of such optimal features defines the dynamics of x p (t) under the general form (which is slightly more complex than (3.4))

dx p dt = l(x p , t)u + s(x p , t) (3.6)
where the nonlinear functions l and s are constructed from the functions T 0 , T px , T p∆ and their derivative or reciprocal function, and parameters S, ρ, c p , T a and k. Simple examples for l and s are reported in (3.4). At any instant t ∈]0, ∆t], the temperature inside the tank is defined as the profile constituted by the plateau (on the lower part) and the initial profile updated by the heat losses factor (on the upper part).

This defines

F II (T (•, t i ), ∆t i , u) = T (•, t i+1 ).

Numerical considerations

In principle, the extra features added to the dynamics could render (3.6) difficult to identify and even more difficult to integrate. However, an integral form similar to the energy balance (3.3) gives an easy way to determine the profile at the end of the heating phase. This method is used in practice to numerically compute the profile in Section 3.3.

Phase IIIa and IIIb: Drain as a Stefan problem

Physical considerations

During draining periods, the model of § 2.1.4 inspired by Zurigat's convection-diffusion model appears to be globally valid at the light of experimental data. However, examination of recordings reveals that the injected water seems to be of higher temperature than the one coming from the water system, and the injection seems to be located not at x = 0 but at higher heights (see Fig. 3.4 (a.2)). Below, we propose an explanation for this.

As we have described it in Section 3.1, the water nozzle mixes the injected water in a neighboring volume, raising the temperature in a zone of varying size. The Zurigat-inspired model does not account for this effect and tends to neglect the water at the bottom of the tank. This results in a undesirable shift of the thermocline (see Fig. 3.4 (a.2)). A similar effect is studied for large storage tanks (having capacity larger than 30 m 3 ) when injecting hot water on top of the tank in [START_REF] Oppel | Computer simulation of stratified heat storage[END_REF] and [START_REF] Nakahara | Initial formation of a thermocline in stratified storage tanks[END_REF]. In these works, buffer zones have been introduced in the proposed models with respectively constant and constantly increasing (with time) lengths. These buffer models do not yield conclusive results in our case, even though interesting similarities in the spirit of derivation can be seen with our work.

For these reasons, we introduce another homogeneous zone characterized with temperature T b (t) and length x b (t). Their dynamics are driven by the water injection. In case of simultaneous heating and draining (case IIIb), draining effects are predominant. We simply assume that the heating elements all belong to this zone and concentrate the effects of the heating power u on T b (t).

Dynamics

We still consider Zurigat's convection-diffusion PDE on the interval [x b (t), h], but with a Dirichlet boundary condition T (x b (t), t) = T b (t) which is now located at the end of the mixing area (x b in Fig. 3.3), and thus constitutes a time-varying boundary condition.

The input variables of F III are the initial profile, say T 0 (•), the duration ∆t i , and the constant values of the heating power u and draining velocity v d . Then, we consider the Stefan problem (3.7)-(3.8)-(3.9)

∂ t T + v d ∂ x T = d α th ∂ xx T -k(T -T a ) on C s T (x b (t), t) = T b (t) on I ∂ x T (h, t) = 0 on I T b (0) = T 0 b , x b (0) = x 0 b T (x, 0) = T 0 (x) on ]x b (0), h]                    (3.7) over the domain C s = {(x, t)|t ∈ I, x b (t) ≤ x ≤ h} (where I =]0, ∆t]). t i Height h x b (t) 0 t i+1 Time T (x, t) = T b (t)
In this domain:

Convection-diffusion equation v d , α d
In this domain: The dynamics of T b derive from the energy balance

x b (t i )
x b (t) dT b dt (t) = v d [T in -T b (t)] + d α th ∂ x T (x b (t), t) -kx b (t)(T b (t) -T a ) + u Sρc p . (3.8)
The dynamics of the moving interface x b (t) are defined as follows. Experimental data suggests that the mixing zone is larger when its temperature is low. Then, we choose the general formulation

x b (t) = q(T b (t)) (3.9)
where q is a positive, continuously decreasing function to be identified. For the simulation, we use the nonlinear expression

q(T ) = min(h, A T -T min ) (3.10)
where A and T min < T in are subject to an identification procedure. More generally, q represents the mixing effects of water injection and therefore is strongly connected with the typology of the nozzle. This definition is consistent with the observation that large values of v d induces rapid decrease of T b and increase of x b . The initial conditions x 0 b and T 0 b are defined as follows. The homogeneous domain is initialized with

x 0 b = argmin{x b |q x b 0 T 0 (x)dx = x b } (3.11) and T 0 b = min{T |q(T ) = x 0 b }. As defined in (3.11), if {x b |q ( x b 0 T 0 (x)dx) = x b } is not empty, x 0
b is unique as the minimum of a closed set (due to the continuity of q and of the integral). An empty set corresponds to the case where very few energy is in the tank and the temperature in the tank is close to T in . We consider then that the mixing zone covers the whole tank and set x 0 b = h. This completes the definition of the dynamics. Remark 1. Generally, existence and uniqueness of such time-varying boundary conditions problem is not straightforward. Here, relating to the so-called Stefan one-phase problem [START_REF] Fasano | General free-boundary problems for the heat equation[END_REF][START_REF] Fasano | General free-boundary problems for the heat equation, ii[END_REF][START_REF] Cannon | The one-dimensional heat equation, volume 23 of Encyclopedia of Mathematics and its applications[END_REF] (which models the temperature of water next to a melting layer of ice and other crystal growth problems), under additional assumptions on q and u, we can prove not only that the problem is well-defined but also that for any time t ∈ I, T (•, t) is an increasing function. This directly follows from the general result of [START_REF] Fasano | General free-boundary problems for the heat equation[END_REF][START_REF] Fasano | General free-boundary problems for the heat equation, ii[END_REF]. Remark 2. If ∃t ∈ [0, ∆t] s.t. x b (t) = h, the previously defined dynamics is halted and is replaced for later times by

h dT b dt (t) = v d [T in -T b (t)] -kh(T b (t) -T a ) + u Sρc p x b (t) = h     
(3.12) defining T (•, t) = T b (t) over the whole tank.

The description above defines the function F III : given a constant draining v d and heating u, F III is the mapping from T 0 (•), ∆t i , v d and u to the solution T (•, s) of (3.7), (3.8), (3.9) at time t i+1 . In other words

F III (T (•, t i ), ∆t i , u, v d ) = T (•, t i+1 ).

Numerical considerations

As in § 3.2.2, this system can be solved numerically with finite difference schemes, jointly with an ODE solver for the state T b with the same time-step.

Summary of the multi-period model

The definitions of functions F I , F II , F III above allow an easy computation of the temperature profile at any final time t f : given a succession of phases separated by times T = (t 0 , t 1 , ..., t m = t f ) and an initial profile T (•, t 0 ), intermediary profiles T (•, t i ) are computed by successive applications of these functions, using the constant values of u and v d during each phase.

Remark 3. Denote E ⊂ L ∞ (Ω) the set of increasing, piecewise continuous and continuously piecewise differentiable functions. Each of F I , F II , F III maps the profiles from E to E, ensuring that our problem is well-defined.

Comparison against experimental data and physics-based model

The multi-period model has been compared against experimental data under the same conditions as in Chapter 2. Comparisons with previously obtained numerical results (using the model from § 2.1.4) show clear improvements both during draining periods (with the vanishing of the temperature and thermocline shifts) (see Fig. 3.4 (a.2)(b.2)) and during heating periods (with an improved estimation of the plateau distribution and the backward flow) (see Fig. 3.4 (a.1)(b.1)). The changes have a positive impact on the accuracy of the model, as is reported in Table 3.1, in which the statistics of the absolute difference between simulation and experimental data has been used as a quality index (for the whole set of data): with the multi-period model, 96.8% of predicted temperatures have an error lower than 4 o C. Most importantly, a significant speedup is obtained. The low computational complexity of the multi-period model enables its implementation in the envisioned "smart EHWT" applications.

Err. 0-2 o C 2-4 o C 4-6 o C 6-8 o C

Summary

In this chapter, we have proposed a new model for an EHWT. This model is based on experimental observations. It has the advantages of being accurate and computationally light.

Thanks to these performance improvements, the multi-period model can be used to address the optimization problems P 1 , P 2 considered in the thesis.

Very generally, in the context of "smart EHWT", this model can be used to estimate the temperature profile of the water contained in the tank, from an initial profile, heating and draining histories. In practice, this requires a flow sensor (either at the input or output of the tank), and means for logging the power injected via the heating element. To estimate the initial condition, an internal temperature sensor can be used. In most EHWT installed in the French market, such a sensor is already installed at the bottom of the tank in the thermostat. It is coupled with the heating element (see Fig. 1.1): heating goes on until this sensor detects a temperature of T max . The description of the thermostat, and the properties of the two models introduced earlier suggest that, when the sensed temperature equals T max , the temperature of the water is uniformly equal to T max .

When coupled with the flow sensor, the thermostat sensor allows us to identify an openloop model with occasional resetting of the profile. Starting with a uniform temperature profile T max , we compute the subsequent profiles using the model and the data of the injected power and drained water. At each occurrence of T max measured by the sensor, we replace the profile by an uniform profile to ensure that the profile is calibrated on a regular basis. These can form the basis of a general state observer for this distributed parameter system.

Chapter 4

Autoregressive modeling for domestic hot water consumption

Modèlisation autoregressive de la consommation d'eau chaude sanitaire. Dans ce chapitre, nous étudions les propriété statistiques de la consommation d'eau chaude sanitaire au cours du temps. Sont étudiés : la distribution des amplitudes des soutirages, le comportement du consommateur sur une journée, et les durées entre deux soutirages successifs. Un résultat marquant est que la durée entre deux soutirages successifs suit une distribution de Weibull. Un modèle conditionnel autorégressif est proposé pour la suite de ces durées, il est validé par des données expérimentales.

In this chapter we study one of the inputs of the previously introduced model: the consumption of domestic hot water.

Users drain hot water from their EHWT at various times of the day. For comfort, delivery of hot water at all times is a constraint that must be satisfied at best. In advanced heating control strategies, a model of the user demand can be used, through the multiperiod model developed in Chapter 3 to estimate, at any time, the state of charge of the tank as a function of heating and domestic hot water consumption. This states serves to evaluate (and guarantee) the availability of hot water for the user. In this chapter, we perform some investigations on the dynamics of domestic hot water (DHW) consumption, which has a random nature.

In the literature, many studies have focused on describing weekly or seasonal consumption patterns, or have presented hour-per-hour mean consumptions [VDS87, MT97, PPS95]. At smaller times scales, the water drains appear as a sequence of quasi-instantaneous drains (see Fig. 4.1). As can be observed on experimental data, the time of occurrence of these drains is not fixed, but is stochastic. Some models representing drains as point process have been developed. One such typical approach can be found in [START_REF] Jordan | Influence of DHW load profile on the fractional energy savings: a case study of a solar combi-system with TRNSYS simulations[END_REF], in which an approach based on aggregation of types of uses allows to generate minute-per-minute load profiles, in [START_REF] Prud | Advanced control strategy of a solar domestic hot water system with a segmented auxiliary heater[END_REF], in which forecast over 2 days are generated using the Kalman filter, and in [START_REF] Hendron | Tool for generating realistic residential hot water event schedules[END_REF] which distribute clusters of drains following a hour-per-hour mean consumption pattern. Nevertheless, more advanced stochastic modeling of the temporal correlation in the sequence, based on data analysis, would represent some valuable ways of improvement, especially in our context of heating strategy design. This chapter aims at establishing such statistical properties. We develop an autoregressive model for domestic hot water consumption. Domestic hot water consumption aggregates various uses in a household: baths, cleaning, cooking, etc. The durations of resulting drains range from a few seconds to a few minutes. They can be represented as quasi-instantaneous events of various magnitudes in the scale of the day. These drains can be described as a volume of hot water or a quantity of energy taken from the tank (the energy contained in hot water being defined with respect to a cold water temperature reference). To model the drain sequence, one needs to describe when the drain happens, and how much hot water is consumed during those drains (we make the assumption that two drains occurring within a single minute correspond to a single bigger drain). The following statistical properties are considered: i) the distribution of the magnitude of the drains, ii) a daily pattern of the start times of consumptions, iii) the time between two successive drains. This defines a stochastic process for the total hot water consumption DHW C(t) at time t of the form

DHW C(t) = +∞ j=1 M j θ(t -t j ) (4.1)
where t j and M j are the time of occurrence and the magnitude of the drain j, resp., and θ(t) is the Heaviside function: θ(t) = 0 if t < 0 and θ(t) = 1 if t ≥ 0. Classically, the magnitude of the drains can be characterized through frequentist inference [START_REF] Neyman | Outline of a theory of statistical estimation based on the classical theory of probability[END_REF], by estimating a probability density function from the frequency of the data. Results obtained with this approach are presented in Section 4.1. The times of occurrence of the drains are more complex to describe. The consumption start times are related to the number of persons in the household and their domestic habits. As will be shown in Section 4.2, they are distributed according to an average daily pattern. Then, the time between two successive drains follows a Weibull distribution. This is shown in Section 4.3. This observation allows to define a Weibull Autoregressive Conditional model for DHW consumption, which is presented in Section 4.4. There, a validation is performed using data gathered in 11 distinct households with EHWT equipped with flow meters and temperature sensors, over periods ranging from 292 to 337 days, with a sampling time of one minute.

Frequentist inference for the water drains distribution law

A simple way to model the diversity of magnitude of the drains is to represent the drains as random variables governed by the same probability law. Since the drains are positive, this probability law can be represented with a probability density function on [0, +∞[ (see Fig. 4.2). The probability distribution reported in Fig. 4.2 has been estimated with a frequentist approach, for all households and all measurements. The displayed distribution law is the frequency of occurrence of each drain for a large number of measurements. An interesting observation is that four distinct peaks are visible in the filtered distribution function. This is in accordance with the results presented in [START_REF] Jordan | Influence of DHW load profile on the fractional energy savings: a case study of a solar combi-system with TRNSYS simulations[END_REF], obtained by considering four types of uses (small and medium drains, shower bath and bath tub). If needed, these results could be refined by taking into account that the measurements are strongly related to the habits of the households (e.g. baths), the time of the day and the season. For households equipped with additional sensors, the conditional distribution could be estimated online. However, we do not consider this refinement here.

Daily pattern for the start times of drains

The hot water consumption is strongly related to the habits of the house occupants. Also, it is clearly visible in the data (as could be expected) that the drains mostly take place at certain times of the day (see Fig. 4.1). In fact, a daily pattern can be defined. For this, one can simply consider the mean value of the number of drains at a given time of the day, for each household. The procedure is as follows. This procedure is commonly used in statistics (e.g. to determine the intensity function of a non-homogeneous Poisson process [START_REF] Ruwet | Processus de Poisson[END_REF]). Given a household, we represent each of its n set of 24 h measurements (labeled by k = 1...n) with a function N k (t) that represents the number of drains having occurred in this day over the interval [0, t], for t ∈ [0, 24 h]. In other words, consider any day k and the sorted drain times t k 1 , ..., t k m , then N k (t) is defined as

N k (t) = m j=1 θ(t -t k j ). (4.2)
Then, we define the daily pattern as the mean of these functions

M (t) = 1 n n k=1 N k (t). (4.3)
This function, defined over the time interval [0, 24 h], gives, for any household, the average (expected) number of drains that should have occurred at time t (see Fig. 4.3): it is representative of the habits of this household in terms of frequency of the drain. The

Distribution of the time between two successive drains

The time between two successive drains y k i = t k i -t k i-1 (during a day k) is also related to the habits of the house occupants. Therefore, we isolate data accordingly. In the case when the number of drains over time N k (t) are samples of a non-homogeneous Poisson process (i.e. the drains are not correlated), it is possible to construct adjusted increments with the mean value function M (t), that must follow an exponential distribution. As will appear, this assumption is not valid in our case, but we still perform the same analysis, yielding different conclusions. The construction is done the following way. For each period k of 24 h, a set is constructed by taking the image of each drain time t k 1 , ..., t k m , through the function M (see Fig. 4.4). The M (t k i ) are then used to define the set of successive increments

υ k i = M (t k i ) -M (t k i-1
). These variables have an interpretation: for any day k, υ k i is the average increment of the number of drains from t k i-1 to t k i . Its expected value (for all possible k) is 1. These normalized increments (represented in the ordinate axis in Fig. 4.4) are assumed to be independent from the choice of the households. They are representative of the correlation between successive drains. The υ k i are distributed in the [0, +∞[ interval. Using the whole dataset, an experimental cumulative distribution function (CDF) can be obtained. It is reported in Fig. 4.5. Remarkably, this function shows a good fit to h(υ) = 1 -exp(υ δ /κ), which corresponds to the CDF of the Weibull distribution of shape parameter δ and scale parameter κ [START_REF] Murthy | Weibull models[END_REF]. Among other possibilities (exponential distribution, Gamma distribution, Beta distribution), this is (by far) the best fit. In the literature, the Weibull distribution is commonly used to model failure rates over time, the parameter δ defining the nature of the process. The case δ = 1 corresponds to an exponential distribution. If δ < 1, then the occurrence of an event raises the probability of a closely following event (cases of immediately consecutive events are frequent). On the contrary, the case δ > 1 corresponds to the case when successive events are spaced out.

To establish that the experimental data follow the assumed distribution, a graphical method can be employed. The probability plot consists in plotting two CDF, one against the other: if the two distributions are similar, the points of the probability plot should lie on a straight line. Further, for the Weibull distribution, the shape parameter δ can simply be deduced from the slope [START_REF] Murthy | Weibull models[END_REF]. Such a plot is reported in Fig. 4.6, where the distribution of the experimental adjusted increment set is plotted in logarithmic scale against a generated data set following a normalized Weibull distribution. The apparition of a straight line in the plot suggests that the Weibull model is relevant. Typically, δ = 0.53. However, the linearity is not perfect for small increments, and a jerky character of the line can be observed for small increments. A reduction of the sampling time can increase the accuracy (here the one minute sampling time is in conflict with the very shorts time intervals separating two successive drains).

Normalized increments

An autoregressive model for domestic hot water consumption

Given the studies above, we incorporate the Weibull distribution in a point process model for the drains. In the literature, the Weibull distribution is often associated with Autoregressive Conditional Duration (ACD) models [START_REF] Crowley | Point process models for multivariate high-frequency irregularly spaced data[END_REF], initially developed by Engle and Russels [START_REF] Engle | Autoregressive conditional duration: A new model for irregularly spaced transaction data[END_REF] 1 . ACD models are based on an autoregressive structure for the conditional duration of time lapses between two events [START_REF] Engle | Autoregressive conditional duration: A new model for irregularly spaced transaction data[END_REF]. The idea is to generate a sequence of n durations with an autoregressive law depending of preceding duration. The duration number i, seen as a stochastic variable υ i , is given by the mean equation

υ i = ψ i i ,
in which the i are independent and identically distributed following a given probability law, with mean value equal to 1, and

ψ i = E[υ i |F i-1 ]
is the mean duration conditionally to the information available after the last event represented by the filtration F i-1 . A filtration is an increasing (in the set theory sense) sequence of σ-fields. It is used to represent the (increasing) information available at a each step. Here, the F i is the natural filtration generated by the sequence of the i . For more information on filtration, see e.g. [START_REF] Rogers | Diffusions, Markov Processes, and Martingales[END_REF].

The linear ACD(1, 1) model is given by the equation

ψ i = ω 0 + ω 1 υ i-1 + ω 2 ψ i-1
for ω 0 > 0, ω 1 , ω 2 ≥ 0. Several distribution can be given to the i . The most classical are the Exponential ACD (EACD) and Weibull ACD (WACD). ACD(p, q) models can be defined by extending the autoregressive character to order p for the υ i and q for the ψ i , but we will limit ourselves to the case (p, q) = (1, 1). For a given household, identification of parameters ω 0 , ω 1 , ω 2 and shape parameter δ of the Weibull distribution is needed. For this problem, we use the off-the-shelf Estimation and Simulation of ACD models Matlab code produced by M. Perlin (see [START_REF] Perlin | Estimation and simulation of ACD models in Matlab[END_REF]).

Then, the identification procedure to model drain occurrence for a given household is the following.

1. Gather historical data on the occurrence time of the drains in the the household.

2. Using this data, compute the daily pattern M (t) and the distribution law of the magnitude of the drains as described in Section 4.1 and Section 4.2 .

3. Using the occurrence time and the daily pattern, construct the set of increments υ i .

4. Identify parameters ω 0 , ω 1 , ω 2 and δ (for instance with the software package [START_REF] Perlin | Estimation and simulation of ACD models in Matlab[END_REF]).

1 It was initially used to model seismic aftershocks and to forecast the distribution of a succession of duration times between two transactions in stock market [START_REF] Tsay | Handbook of Econometrics II, chapter Autoregressive Conditional Duration Models[END_REF] This procedure permits to represent each household with model parameters

M (•), ω 0 , ω 1 , ω 2 , δ.
The WACD parameters are used to generate sets of increments υ i , which are then used to construct time of occurrence of drains using the relation

t i = M -1 (υ i + M (t i-1 )).
Finally, the magnitude of the drain can be drawn from the distribution law previously defined.

Model validation

To validate the model, we apply the identification procedure for each households of the dataset described earlier. Then, the obtained parameters are used to generate sets of increments. Finally, both experimental and simulated increments are compared in a probability plot in Fig. 4.7. The results are satisfactory, at the exception of long experimental durations that are not well represented in the model. A possible explanation is that those duration may correspond to vacation periods (holidays), and represent outliers that we could not isolate from the dataset. This stresses the necessity to add an overlay to account for exceptional events in future studies. As a mean of comparison, probability plot against duration stemming from uniform distribution of drains on a 24 h period, and against an EACD model are depicted in Fig. 4.8 and Fig. 4.9. The conclusion of this analysis is that the Weibull distribution is, by far, the best standard distribution for an ACD model at the light of the data under consideration.

Summary

In this chapter, an autoregressive time series model for DWC consumption has been developed. Interestingly, the model stresses that if the user is draining water, then he is very likely to drain more water shortly. Consequently, an online computing of the conditional mean duration ψ i gives a representative index of the drain general trend. A small value indicates a general increase of consumption, and piloting strategies can be updated online to take into account the high probability of future drains. Very generally, this model can be used for "smart EHWT" piloting applications. Indeed, it can be used to generate realistic scenarios of consumption, subsequently used to test the robustness of piloting strategies. At large scales, for a large number of EHWT, stochastic modelling of the consumption can be generalized into models for the distribution of the state of charge of the tanks, using e.g. Fokker-Planck equations. This will be shown in Chapter 8. To design control and optimization advanced strategies, we define key indicators representing the state of each EHWT. By proceeding this way, we reduce a distributed state (the temperature profile in the tank) into a finite number of well-chosen variables directly useable to evaluate comfort constraints. A focus is put on the supply of water to consumers when piloting the heating of the tank. Considering that the temperature profile defined above is increasing with height because of stratification (see Chapter 2 and 3), further simplifications are possible. One can define two variables representing the energy contained in the water of temperature above a comfort temperature set by the user (i.e. an energy which can be used for consumption), and a energy that represents the functioning delay of EHWT for hot water production. The delay energy represents the energy necessary for the plateau introduced earlier to reach the comfort temperature. Once this level is reached (after the so-called "delay"), further heating raises the energy available for consumption. A third variable of interest is used to represent the energy in the tank unavailable for consumption, i.e. contained in water heated below the comfort temperature.

This triplet of variables offers a simple representation of an EHWT. Interestingly, the complex variation of the temperature profile is well represented by the dynamics of this low-dimensional state. This representation is an intermediate solution between the distributed parameter models developed earlier in the thesis and single integrator models commonly found in the literature. It models the complex dynamics of the comfort variables and the lag observed in the production of hot water, in a simple yet non-simplistic manner.

Definition of user comfort

Following the description above, the temperature T of the tank is a continuously increasing function of height. We assume that the water injected in the tank is at constant temperature T in , which constitutes a lower bound of the temperature profile, and that the heating process is driven by turbulence generated by buoyancy effects during the heating process, which is the cause of a local mixing in the bottom of the tank. We consider that this mixing is perfect on the plateau, and does not affect the temperature profile in the upper part of the tank. During the process, the plateau grows and gradually covers the whole tank. Thanks to a thermostat, the user specifies a temperature T max at which the heating has to be stopped to prevent overheating. As a result of the heating process, if the temperature at the bottom of the tank is T max , then the temperature in the tank is uniformly at T max . The user also specifies a comfort temperature T com . Water having temperature higher than T com can be blended with cold water to reach T com and is therefore useful, while water having temperature lower than T com is useless.

Definition of variables of interest: available, delay and reserve energies

We simplify the system using the 3 following state variables, depicted in Fig. 5.1, The available energy a is defined as the energy contained in the zones having temperature greater than the comfort temperature T com . This constitutes a direct comfort index for the user. The situation when a reaches the value 0 and a water drain is applied corresponds to the fact that the consumer is trying to consume hot water when none is available, and therefore that the comfort constraints are broken.

The delay energy τ is defined as the energy required by the plateau to reach the temperature T com . If the tank is heated at constant maximum power, and without accounting for drains and heat losses, τ is proportional to the time necessary to reach a state in which a can effectively be increased by the heating process.

The reserve energy µ is defined as the energy contained in the tank that is currently unavailable for consumption, i.e. the energy contained in the water under T com . When, during the heating process, τ reaches the value 0, the energy µ becomes available to consumption: this generates an immediate (discontinuous) increase of a, and µ is reset to 0.

The rationale behind these definitions is that to plan the heating, we account for the time left before the energy reserve embodied by a (in the total energy a + µ) is consumed, and the time necessary to provide new hot water, embodied by τ . Fig. 5.2 reports examples of the dynamics of these variables. A drain (shown in Fig. 5.2 (a) and (e)(f)) is mainly characterized by a decrease of a and an increase of τ , with a slight raise of µ due to an energy transfer from a. On the other hand, in the heating reported in Fig. 5.2(b)(c)(d), τ decreases at the same rate as µ rises, and after a gap at τ = 0 characterized by an energy transfer from µ to a, a progressively increases.

Given T (•) the (non necessarily strictly) increasing temperature profile of the water defined on [0, h], the above definitions yield the following expressions of a, τ , µ as In the sequence, we will use the notation z (a, τ, µ).

The triplet z is a low-dimensional representation of the tank. This representation comes at the cost of a loss of information in comparison to a distributed profile. However, due to the increasing nature of the profile and the fact that the thermocline is often steep, the triplet gives a satisfactory schematic view of the distributed profile (see Fig. 5.3). Moreover, it captures well the information necessary to design control strategies.

An example of the variation of z over a 24 h run with given drains and heating is presented in Fig. 5.4. In this figure, the variations of the three variables are computed with experimental data (in blue) and compared with simulation with the multi-period model (in red) and a simple single-zone model (in black). The figure highlights the accuracy of the multi-period model that is able to kindly reproduce the evolution of z.

Summary: towards the "smart EHWT"

The definition of the variable of interest concludes the modeling efforts conducted so far. At this stage we can now formulate an input-output representation of the EHWT, presented in Fig. 1.4. These works allows to design the architecture of a "smart EHWT" equipped with a chip set as depicted in Fig. 5.5.

In such tanks, histories of electricity consumption and water drains allow to compute, thanks to the multi-period model, the distributed temperature profile of the water in the tank. Variables of interest derived from the profile are transmitted to an external decision center, which, in return, transmits a heating strategy back to the "smart EHWT". If necessary, this strategy can even be test-proven by the "smart EHWT" itself, using the domestic hot water consumption model and the multi-period model numerically treated by the onboard chipset. Eventually, the strategy can be applied in real-time (up to communication and computation lags).

Remark 4. Alternative architectures are possible. For instance, an intermediate decision center can also be integrated to the chipset to generate heating strategies for the EHWT with respect to a price signal.

If such smart tanks are to be developed and deployed, their use in coordinated strategies have to be designed to solve P 1 and P 2 . The second part of this thesis focus on such problems. Reserve energy (J)

Part II

Control of groups of EHWT Chapter 6

Small-scale groups: discrete-time optimization

Optimisation discrete pour un parc de petite taille. Dans ce chapitre, nous proposons un modèle discret d'évolution de l'état du ballon représenté par les trois variables d'intérêt définies précédemment. Ce modèle permet de formuler un problème d'optimisation du fonctionnement du ballon sous contraintes de confort. Ce problème peut lui-même se ré-écrire comme un programme linéaire en nombre entiers, de telle sorte que la résolution est réalisable avec des logiciels standards. Des résultats numériques sont présentés et comparés avec la stratégie historique d'heures pleines/heures creuses.

In this chapter, we propose an optimal heating strategy for a first use-case focused on a small group of "smart EHWT" (from 1 to 4, typically). Using the variables of interest introduced in Chapter 5, we set up a discrete-time model for the dynamics of the triplet, as well as an optimization framework to solve problems P 1 and P 2 . This framework takes the form of a mixed integer linear program (MILP) for problem P 1 , or a mixed integer quadratic program (MIQP) for problem P 2 . In the case of P 1 , the computation of the optimal solution is light and may be directly embedded in the chipset of a "smart EHWT"(which could embed mildly powerful processor such as ARM Cortex-A7 MPCore).

The chapter is organized as follows. Section 6.1 is dedicated to the discrete-time model. The dynamics serve to formulate P 1 and P 2 in discrete time. Section 6.2 shows how to recast such problems as mixed integer linear/quadratic programs, which yields easy numerical solution. Numerical results are presented in Section 6.3.

Discrete-time dynamics and optimization problems

Discretization and notations

For a given tank, we consider the variation of z = (a, τ, µ) over a finite time horizon that we discretize into uniform time-steps [0, ..., n]. We note (a 0 , ..., a n ), (τ 0 , ...τ n ) and (µ 0 , ..., µ n ) the values of a, τ , µ at each of these time-steps. At each time t ∈ [0, ..., n -1], energy is consumed by the user (through draining) by an amount d t that we consider as known. Energy is introduced via the heating element by an amount u t ∈ [0, u max ]. We divide the later into two parts v t and w t , representing the share introduced in a t and in µ t , respectively. Finally, we define the variable φ t which represents a flow of energy from µ t to a t+1 , taking place in conditions that will be described below.

Dynamics and constraints

Balance equations

For any index t, the dynamics of z t is given by the energy balances

a t+1 = (1 -p)a t -αd t + v t + φ t (A t ) τ t+1 = τ t + pµ t + βd t -w t (B t ) µ t+1 = (1 -p)µ t -(1 -α)d t + w t -φ t . (C t )
The various phenomena described earlier appear in the righthand side of these equations. Heat losses, modeled with an exponential decay, are characterized by the decrease of a and µ at a rate p, the energy from the later contributing to a raise of τ (see Fig. 5.1). The energy consumed by the user during one time step is split between a and µ with a ratio α, 1 -α and affects τ with a coefficient β.

Definition of the heat source terms v t , w t

As has been seen before, the energy is injected at the bottom of the tank via the heating element. As a consequence, the heating has no impact on the available energy a when τ > 0, but, instead, tends to reduce τ and increase µ. When the value of τ is 0, the injected energy becomes immediately available. This can be modeled by dividing u t into two shares v t and w t representing, respectively, the part of the injected energy going into a and µ, and subject to the following conditions

u t = v t + w t (D t ) 0 = v t τ t+1 (E t ) 0 ≤ v t , w t ≤ u max . (F t )
Given the balance equations (D t )-(F t ), if τ t+1 > 0 then no energy can be introduced in a t+1 (i.e. v t = 0 and w t = u t ) and, if τ t+1 = 0 the value of w t has to compensate for heat losses pµ t and/or energy drain from the user βd t , while the remainder is introduced in a t+1 .

Definition of internal energy flow φ t

The flow φ is always equal to 0, except when τ reaches the value 0. Then, the value of φ is defined by the fact that all the energy µ suddenly becomes available. This can be described as follows

0 = φ t τ t+1 (G t ) 0 ≤ φ t (H t ) τ t+1 = 0 ⇒ µ t+1 = 0. (I t )
Remark 5. Given a state z t , a drain d t and a heat injection u t , relations (A t )-(I t ) uniquely define v t , w t , φ t and therefore the future state z t+1 (under the assumption that τ t+1 has to be nonnegative).

Bounds and comfort constraints

Here, we determine some constraints z t is subject to for all t. This defines the admissible controls that allow each z t to respect these constraints. By definition, the delay energy is positive. We impose τ t ≥ 0 for all t, in order to impose a heat injection with v t if τ reaches 0.

Define e max = Sρc p (T max -T in ) the maximal energy that can be contained in the tank, and

λ T com -T in T max -T in . (6.1)
To ensure that no energy is drained more that the tank can provide, we require a t ≥ 0 and µ t ≥ 0 for all t, and that the tank is not overheated, we require a t ≤ e max , τ t ≤ λe max , and µ t ≤ λe max for all t.

Finally, physical constraints on the total energy imply

λa t + τ t + µ t ≤ λe max λe max ≤ a t + τ t + µ t . (6.2)
Given these relations, we define the following domain

Ω = {(a, τ, µ) ∈ R 3 + |λe max ≤ a + τ + µ and λa + τ + µ ≤ λe max }.
Then, ∀t, z t is subject to the constraint

z t ∈ Ω. (J t )

Admissible controls

For any given (a 0 , τ 0 , µ 0 ) ∈ Ω and d = (d 0 , ..., d n-1 ) ∈ R n + , and for any chosen control sequence u = (u 0 , ..., u n-1 ) the relations (A t )-(I t ) for t ∈ [0, ..., n -1] uniquely define (a, τ, µ, v, w, φ) t for all subsequent t. This allows us to define the admissible set U:

U(a 0 , τ 0 , µ 0 , d) = {u ∈ R n |the a, τ, µ, v, w, φ implicitly defined by (A t )-(I t ) t∈[0,...,n-1] satisfy (J t ) ∀t ∈ [0, ..., n]}.
In a practical sense, given initial conditions and a drain sequence, it constitutes the set of heating sequences that do not break the comfort constraints of the user.

Objective function and formulation of the optimization problem

We consider that a 0 , τ 0 , µ 0 and d are given and that a flat discrete-time price signal for electricity over time (c 0 , ..., c n-1 ) is known. The optimization problem P 1 can be written (for each tank) as min

u∈U (a 0 ,τ 0 ,µ 0 ,d) n-1 t=0 c t u t . (6.3)
This problem is therefore to search, for an EHWT, the heating strategy that minimizes the cost of heating while ensuring a required supply of hot water to the user. In cases when several tanks can simultaneously be controlled with controls u 1 , ..., u k , and one can introduce the discrete-time objective for aggregated consumption (P t ) t∈[0,...,n-1] . Then, the optimization problem P 2 takes the form min

(u 1 ,...,u k )∈U (a 1 0 ,τ 1 0 ,µ 1 0 ,d 1 )×...×U (a k 0 ,τ k 0 ,µ k 0 ,d 1 ) n-1 t=0 (P t - k j=1 u j t ) 2 .
(6.4)

Mixed integer representation of the constraints and dynamics

With the exception of the two product conditions (E t ) and (G t ), and the condition (I t ), the relations (A t )-(J t ) are linear equalities and inequalities in the variables (u t , a t , τ t , µ t , v t , w t , φ t ) and therefore define a polytope of R 7n , for each tank. Thus, we adapt (E t ), (G t ) and (I t ) to give to problems (6.3) and (6.4) the structure of a linear and a quadratic program, respectively, that can efficiently be solved with commercial software (see e.g. [START_REF]Gurobi optimizer reference manual[END_REF] or [START_REF] Ilog | User's Manual for CPLEX[END_REF]).

Strengthening of (I t )

Given the set

A = {(τ, µ) ∈ [0, λe max ] 2 |τ = 0 ⇒ µ = 0}, then, equivalently, A = [0, λe max ] 2 \{(τ, µ) ∈ [0, λm] 2 |τ = 0, µ > 0}.
A possible strengthening is

{(τ, µ) ∈ [0, λe max ] 2 |M τ t ≥ µ t } ⊂ A
where M > 0 (see Fig.

Y(a 0 , τ 0 , µ 0 , d) = {Y = (u, a, τ, µ, v, w, φ) ∈ R 7n | (A t )-(D t ), (F t ), (H t ), (I t ) ∀t ∈ [0, ..., n -1]
and (J t ) ∀t ∈ [0, ..., n] are satisfied}.

Reformulation of (E t ) and (G t )

The product conditions (E t ) and (G t ) correspond to 2n couples (v t , τ t+1 ) and (v t , φ t ) consisting of two elements that cannot simultaneously be positive. These constraints have the form of complementarity conditions (see [START_REF] Cottle | Complementary pivot theory of mathematical programming[END_REF]). In practical terms, each of these cases can be encoded with binary or integer constraints and problems (6.3) and (6.4) written as a collection of linear programs (LP), respectively quadratic programs (QP), taking the form of a mixed-inter linear program (MILP) and a mixed integer quadratic program (MIQP).

A more careful look at the situation reveals that these kinds of sets are in fact frequently encountered in discrete optimization under the name of Special Ordered Set (SOS), and are associated with branch and bound strategies (see [START_REF] Beale | Special facilities in a general mathematical programming system for non-convex problems using ordered sets of variables[END_REF]). This considerably eases the resolution for instances having reasonable dimensions. These strategies are often implemented in linear programming solvers, see [START_REF]SCIP User's Manual[END_REF].

As a result of the reformulation of (E t ) and (G t ), and the strengthening of (I t ), we can propose the following solution method.

Solution method

Given initial conditions (a 0 , τ 0 , µ 0 ) and the drain sequence d, one can now solve by any appropriate numerical method (e.g. branching strategies) the collection of LP/QP min

Y ∈Y(a 0 ,τ 0 ,µ 0 ,d) n-1 t=0 c t u t (6.5) or min Y 1 ,...,Y k ∈Y(a 0 ,τ 0 ,µ 0 ,d) n-1 t=0 (P t - k j=1 u j t ) 2 . (6.6)
indexed by the 2nk SOS on τ and µ.

Simulation results

Identification of the parameters

In the previously defined problems, values for most of the parameters are easy to determine: T com , T in , T max and therefore λ are either chosen by the user, or directly measured from the water distribution network. On the other hand, the values of e max , p and u max depend on the type of EHWT. These data are provided by the manufacturer, and are reported in Table 2.1.

The main difficulty is to determine the coefficients (α, β) of our model. They represent the effects of the energy drains (d 0 , ..., d n-1 ) on z.

In theory, the values of these two parameters depend on x, which breaks the linearity. In practice, their variations are small, except in some parts of the domain Ω, where the temperature in a is close to T com (close to the border a + τ + µ = λm). Therefore, to keep the linearity of the constraints, we assume that the two parameters are constant. To reduce the cost of this assumption in terms of accuracy, we replace the constraint (6.2) with a + τ + µ ≥ λe max + e p , with e p appropriately chosen (for instance e p = 1 10 e max ). Using the model presented in Chapter 3 for identification purposes, we set the values α 1.2 and β 0.4.

Simulations results

To test the relevance of our predictions, simulations with realistic parameters have been conducted on simple cases that can readily be interpreted. The parameters correspond to the tank used for experimentation. The computations have been performed using the SCIP solver (see [START_REF]SCIP User's Manual[END_REF]).

Several cases have been considered. We now detail the obtained results.

Two prices signal

In Fig. 6.2, we report the values of the control and the associated available water for problem (6.5) with those prices. The optimization has been performed over a 48 h horizon with 20 min time step for a 24 h application of the control signal, to avoid side effects. The results show a reasonably smooth control signal, which heats the tank at the end of low prices period, to limit consumption at high prices and reduce heat losses. This strategy is more efficient than the heating at the beginning of the off-peak period, usually set up in practice (the so-called night time switch). For this typical case, the gain for the consumer is of 13.2%, due to a heating strategy that limits heat loss while heating at the end of low-price periods, and avoids overheating of the tank.

Spot prices

In Fig. 6.3, we report the values of the control and the associated amount of available water for problem (6.6) with the spot prices in Europe for a typical winter day, taken from the EPEX Spot Market (see [START_REF]Annual report[END_REF]). Again, the results show a sensible control, which heats the tank at the end of low prices period, to limit consumption at high prices and reduce heat losses. The gain for the producer, in this case, is of 36.1% compared to the night time switch strategy.

Peak shaving for several tanks

As an extension, we can also solve a problem of practical interest. We have run simulations on a set of several tanks with various initial conditions, with a problem defined as weighted sum of the criteria of problem (6.5) and problem (6.6): each tank is optimized against a price signal, individually, and a global consumption is favored in the middle of the day in the form of an additional load curve objective in the [10 h, 13 h] time window that may, for instance, answer to an upward adjustment request. The results of these simulations are reported in Fig. 6.4, where each tank is represented with a different color.

Summary

This chapter has proposed numerical methods to address the optimization problem P 1 and P 2 at stake in the thesis. The method is capable of handling a single or a small group of EHWT. For this a simple dynamics for the 3 states representation of the EHWT, and a framework to optimally heat it with off-the-shelf solvers has been derived. We have processed the optimization problems in deterministic cases, under perfect knowledge of consumptions over the optimization horizons. Practical online computation of heating strategies should account for the inherent stochasticity of users consumptions. The definition of a set of realistic scenarios of drains for a given household using the model presented in Chapter 4 allows one to define a worst-case heating strategy. Stochastic programming also appear as promising solution but has not been explored.

If necessary, some sort of closed-loop feedback model predictive control could be employed. Typically, in the envisioned "smart EHWT", heating strategy could be computed on a 48 h horizon, with 15 min closed-loop intervals.

The method scales reasonably up to a number of 4 tanks. Typical computational loads are reported in Table 6.1. Generalization to arbitrary numbers of tanks does not seem feasible, as-is. The computational cost seems to increase exponentially with the number of tanks. Thus, a global piloting using a single centralized optimization program does Number of tanks Computational time (s) 1 2.0 2 9.4 3 754.2 4 4756.4

Table 6.1: Computational times.

not appear as a viable solution. However, the tanks having separable constraints, the optimization problem under consideration is compliant with decomposition-coordination techniques. In particular, the so-called price coordination (see [START_REF] Mesarovic | Theory of hierarchical multilevel systems[END_REF]) seems to apply well to the formulations considered in this chapter.

Chapter 7

Medium-scale groups: optimization algorithm for load-shifting Algorithme d'optimisation du décalage de charge pour les parc de taille intermédiaire. Dans ce chapitre, nous proposons la formulation du problème P 2 dans le cas ou les périodes de chauffe restent indivisées, tout en prenant en compte l'effet des pertes thermiques. Dans un second temps, une heuristique pour ce problème est proposée. Cette heuristique est testée sur huit instances représentatives des problèmes rencontrés en pratique.

In this chapter, we focus on a medium-sized groups of EHWT (from a few hundreds to several millions), each being characterized by its electric power and heat losses coefficient. We assume that they are remotely controllable.

As described in the introduction of the thesis, a consequence of the night time switch using wired remote control of the group of EHWT in France is an overall decrease of the consumption in the middle of the night when the electricity production costs are low. This is a detrimental effect for the electricity producer.

Instead of finely piloting each tank, we consider a rescheduling of each EHWT heating in the night time, while keeping the heating period undivided, to attain a desired load curve for the whole set of EHWT, while ensuring individual users comfort. The idea is to minimize the malicious impact of complex dynamic effects taking place inside each EHWT described in Chapters 2 and 3, while heating during periods of low general consumption, and to propose a robust strategy which does not divide the heating during the day.

The first contribution of this chapter is to formulate the rescheduling of the set of EHWT heating with comfort constraints as an optimal control problem of the form P 2 , in the case where every tank is remotely controllable and smart (i.e. capable of computing and transmitting information on their current state to the decision center). The optimization problem considers hot water consumptions, production objectives and comfort constraints. Initially, a certain heating starting time is assigned to each EHWT having a uniquely defined duration for the subsequent indivisible heating period. Each duration is scaled according to next day hot water consumptions. The sum of all the power consumptions defines the initial load curve, as a function of time. The optimal control formulation is the first contribution of the chapter. It is discussed in Section 7.1.

A heuristic specifically designed to solve the considered discrete-time quadratic for-mulation is proposed. This is the second contribution of the chapter. The design of any such heuristic strongly depends on the distribution of the durations of the un-dividable heating. Here, a reference distribution obtained from data measured in a vast set of French households is presented and studied. The dataset contains heating periods lasting up to 8 h. The definition of the proposed heuristic stems from the following considerations. First, one notes that the scheduling of EHWT with long heating durations is the most likely to generate high consumption in on-peak period. Secondly, a rigid individual scheduling for a high number of EHWT is prone to generating singularities resulting in undesirable high consumption peaks in the load curve. The heuristic is designed to reduce singularities while generating a high diversity in the distribution. For this purpose, we introduce some stochasticity. Its governing principle is to sort the EHWT by decreasing order of duration times, and, then, to randomly schedule them one-by-one according to an adaptive distribution law. Each duration is compensated to account for heat losses. For a given EHWT, this distribution law depends on the residual load curve, which is obtained by constructing the objective load curve minus the power consumption of the tanks already scheduled. To maximize diversity, the distribution law favors scheduling within time periods containing only few previously rescheduled EHWT. Optimization results produced with real data are given. Objective load curves have been provided by the French utility EDF, and several distribution laws have been tested. They lead to an optimality loss of less than 1%. In a second phase, we extend the formulation of the optimization problem in the case when only a fraction of the set of EHWT is smart. For these, noisy (unbiased) estimates are assumed to replace the missing information. The formulation intends to make use of the flexibility of each tank while ensuring comfort even for tanks from which few information is available from. The same heuristic is used to produce a second round of simulations, in which the effects of uncertainty are quantified.

The chapter is organized as follows. In Section 7.1, the optimization problem is formulated. Then, in Section 7.2 , the heuristic is presented. Simulation results are reported in Section 7.3, while Section 7.4 focuses on the extension to the uncertain case. Conclusions and perspectives are given in Section 7.5.

Formulation of the problem

Nature of the problem

In this section, we consider a pool of n tanks. Each tank starts with a given energy content, and has to be heated during an undivided period of time to reach a given final energy content (which is the one it will have reached with the night time switch strategy). The energy in each tank is subjected to heat loss.

Here, we formulate a problem taking the form of P 2 that includes these requirements.

Electric water heating

We represent each tank (labeled with i ∈ [1, ..., n]) using the energy e i (t) = a i (t) + µ i (t) of the water contained in the tank at time t (defined relative to the energy of cold water). This energy ranges between 0 and a maximum value e i m . It shall be noticed that e i m depends on the volume of the tank, and also on the maximal acceptable temperature in the tank T max set by the user. We lump energy losses to the ambient into a heat losses coefficient k i , Hour representing the heat losses per unit of time (as a percentage of the total energy). We note u i max the power of the heating element. We assume that this power cannot be modulated, i.e. that if the tank is heating, then it is with the (constant) injected power u i max . We focus on a time interval [t 0 , t f ] in which the tanks are heated, usually during the night. For each tank i, we define a time interval [t i 0 , t i f ] ⊂ [t 0 , t f ] during which heating the tank i is allowed. For instance, in France, households benefit in the night-time of block of hours during which the electricity price is reduced to promote electricity consumption when the prices are low, and EHWT heating has to be performed only during these off-peak periods. Further, we define for each tank i the time t i c of the first (earliest) consumption of hot water, at which a certain quantity of hot water must be contained in the tank. The state-of-the-art heating policy (so-called night time switch) is to switch-on each EHWT at the beginning of its off-peak period until it is fully heated. This strategy leads to an aggregated load curve that rapidly decreases in the middle of the night, when the prices of electricity are the lowest. Depending on the season and the market prices, various load curves can bring substantial savings for the electricity producer (see Fig. 7.1 where peak consumption of EHWT is represented along market prices during a typical day). The application of the night time switch does not necessarily generate such a desirable load curve.

Generally, one can describe one of the goals of the electricity provider as follows: given a reference load curve f a : [t 0 , t f ] → R + corresponding to the heating of the tanks from the energy e i (t 0 ) = e i 0 to e i (t f ) = e i f for all i describing the group of tanks, how can all the heatings be rescheduled to approach an objective load curve

f o : [t 0 , t f ] → R + ?
Total power Time

t 0 t f t 0 + ∆t a t 0 + ∆t b d a d b < d a

Preliminary: time of heating compensation and feasibility of load curves

Any rescheduling of the heating of any EHWT has to take into account comfort constraints. Due to heat losses, the duration of heating must compensate changes in the starting time.

If the EHWT is heated sooner (respectively later), heat losses are increased and more (respectively less) energy has to be injected into the EHWT. For any EHWT, a change in heating time can be analytically computed as a function of its heating starting time. We address the effect in the following way.

For any given tank i ∈ [1, ..., n] that, in the so-called reference scenario, starts its heating at time t 0 + ∆t i a for a period of length d i a , we wish to estimate the heating duration d i (∆t i ) that is needed if one chooses to start the heating at time t 0 + ∆t i , without altering the energy of the tank obtained at final time t f (see Fig. 7.2).

If we consider that the ambient temperature T a is equal to the injection temperature T in , space integration of the temperature in the tank given by the dynamics in Section 2.1 or 3.2 leads to de i dt = -k i e i + u i -c i , where c i (t) is the energy consumption in the tank at time t (assumed to be independent from the state). Then, Duhamel's formula yields

e i f = exp(-k i (t f -t 0 ))e i 0 + t f t 0 exp(-k i (t f -s))u i max 1 [t 0 +∆t i a ,t 0 +∆t i a +d i a ] (s)ds -C i (7.1)
and

e i f = exp(-k i (t f -t 0 ))e i 0 + t f t 0 exp(-k i (t f -s))u i max 1 [t 0 +∆t i ,t 0 +∆t i +d i (∆t i )] (s)ds -C i (7.2)
where 1 [x,y] (s) = 1 if s ∈ [x, y] and 0 otherwise, and

C i = t f t 0 exp(-k i (t f -s))c i (t)ds
is the total energy consumption on [t 0 , t f ] for the tank i. Subtracting (7.2) to (7.1) leads to

∆t i a +d i a ∆t i a exp(k i s)ds = ∆t i +d i (∆t i ) ∆t i
exp(k i s)ds, (7.3) and, finally,

d i (∆t i ) = d i a + 1 k i ln(exp(k i (∆t i -d i a )) + exp(k i ∆t i a ) -exp(k i (∆t i a -d i a ))) -∆t i . (7.4)
Then, if ∆t i > ∆t i a , we have d i (∆t i ) < d i a , and if ∆t i < ∆t i a , we have d i (∆t i ) > d i a (see Fig. 7.2).

This delay effect has an impact on the global load curve. Indeed, to reschedule the time of heating of tanks from a reference load curve f a to an objective load curve f o , one needs to ensure that the energy injection is compensated. For a set of tanks with initial starting times t 0 + ∆t i a and heating durations d i a , we have

f a = n i=1 u i max 1 [t 0 +∆t i a ,t 0 +∆t i a +d i a ] .
If we consider that the heat losses coefficients of the tanks are all close to their average k (i.e k i k), the same reasoning as before leads to

n i=1 e i f = exp(-k(t f -t 0 )) n i=1 e i 0 + t f t 0 exp(k(t f -s))f w (s)ds - n i=1 C i
for w = a, o and then to the following "feasibility" condition for f o t f t 0 exp(ks)(f o (s) -f a (s))ds = 0. (7.5)

Problem formulation

In this first case, we assume that all the tanks are smart. In this case, for each tank i the values of e i (t 0 ) = a i (t 0 ) + µ i (t 0 ), and of τ i (t 0 ) are known. To guarantee the availability of hot water during the earliest consumption, we consider the additional constraint that, at this time, the plateau formed during the heating process must have already reached the comfort temperature. This constraint, for a heating starting time ∆t, takes the form ∆t + exp(k i ∆t)τ i (t 0 ) ≤ t c because at rest, the delay energy increases exponentially at rate k i as has been seen in § 6.1.2. We now define our optimization problem. Given n EHWT characterized by their heat losses coefficients k i , their powers u i max , their allowed time intervals [t i 0 , t i f ], their first consumption time t i c , and given an initial load curve f a and an objective load curve f o verifying (7.5), we desire to solve

min ∆t 1 b ,...,∆t n b t f t 0 ( n i=1 u i max 1 [t 0 +∆t i b ,t 0 +∆t i b +d i b ] (s) -f o (s)) 2 ds (7.6) s.t. ∀i = 1, ..., n, the couple (∆t i b , d i b ) satisfies d i b = d i (∆t i b ), (7.7)
with d i as defined in (7.4),

t 0 + ∆t i b ≥ t i 0 , t 0 + ∆t i b + d i b ≤ t i f , (7.8) and ∆t i b + exp(k i ∆t i b )τ i (t 0 ) ≤ t c .
(7.9)

Solution method

We propose here a resolution method for (7.6) under the constraints (7.7)-(7.8)-(7.9), with a discretization of the time window [t 0 , t f ] into p time-steps. In discrete time, the problem formulation presented in § 7.1.4 does not change, but starting times have to be chosen among the p possibilities. Mathematically, the problem defined above can be related to capacity scheduling problems [START_REF] Veinott | Optimal capacity scheduling[END_REF] or cumulative non preemptive scheduling [START_REF] Carlier | Jackson's pseudo-preemptive schedule and cumulative scheduling problems[END_REF].

However the large number of optimization variables, the flexibility yielded by the cumulative nature of power1 are not accounted for in classical resolution techniques. Moreover, the problem is also specific: i) the power consumed in one tank is very small compared to the whole set, because the set is large; ii) the duration of the longest heating (up to 7 h) is relatively bulky compared to the time horizon (10 to 12 h). In discrete-time, determining whether such a problem admits an optimum solution is equivalent to the classical exact cover problem [START_REF] Jongen | Optimization theory[END_REF], which consists in exactly covering a set with a sub-collection of its subsets. Due to the high complexity of solving this problem (shown to be NP-complete in 1972 by R. Karp [START_REF] Karp | Reducibility among combinatorial problems[END_REF]), and more generally, to the difficulty of minimizing such a criteria with a large number of decision variables (the number of tanks can be up to several hundreds of thousands, or millions), we propose in this section a heuristic in the form of a stochastic sequence of rescheduling.

The proposed heuristic is based on the following observations:

• The tanks with long heating durations are more difficult to allocate. On the contrary, short durations yield flexibility to our problem.

• A rigid rescheduling of heating times according to a deterministic procedure is prone to generate singularities in the resulting load curves. Indeed, depending on the shape of the objective load curve, a deterministic rescheduling is very likely to create patterns which generate undesirable high consumptions peaks. On the contrary, due to the large number of tanks, a stochastic heuristic can take advantage of the natural smoothing induced by introducing diversity in the reallocation.

The heuristic we propose is based on successive updates of a residual load curve f i r (t) representing the objective load curve minus all powers from the tanks that have been already rescheduled (see Fig. 7.3). The steps are:

1. Compensate all durations as if the tanks were all starting heating at t 0 , using equation (7.4).

2. Sort the tanks by decreasing compensated durations.

3. Initialize f 0 r (t) = f o (t). 4. For all tanks, from i = 1 (which heats the longest) to i = n (which heats the shortest), apply the following steps (a) Using f i-1 r (t), define a (finite) set of admissible starting times S i , which is the set in which the starting times can be chosen. In practice, S i ⊂ [t i 0 , t i f -d i (0)], and overloaded periods are excluded from the set. (b) Using f i-1 r (t), define a (discrete) probability distribution law L i on S i that promotes rescheduling in under-loaded periods.

L i . (d) Update f i r (t) = f i-1 r (t) -u i max 1 [t 0 +∆t i b ,t 0 +∆t i b +d i b ] (t), s.t. d i b = d i (∆t i b ).
In the heuristic above, the choices of S i and L i are important. In practice, we give S i the form of a union of disjoint intervals (pictured in green in Fig. 7.3 (a) and (b)), corresponding to starting times ∆t i for which the duration d i (∆t i ) is entirely included in times such that f i-1 r (t) > 0 (see Fig. 7.3 (b)). This allows to prevent scheduling in already fully-loaded periods. Several types of laws can be considered, and have various efficiencies depending on the shape of the objective load curve. From simulations performed to investigate this point, we note that if very high slopes are present in the objective load curve, heavy weight should be placed on the boundary of the interval, to promote rescheduling near the boundaries. For this purpose, a parabolic distribution law can be proposed. The law can also directly be defined as a weighted integral of f i-1 r (t) (see Fig. 7.3 (c) and examples of distribution laws in Fig. 7.4).

Simulations results

Dataset

The efficiency of the heuristic strongly depends on the shape of the objective curve, the time intervals, and the flexibility yielded by the diversity of the durations among the population of the tanks. For this reason, realistic data have been gathered for testing. The distribution of compensated durations has been constructed based on data gathered on a panel of 267 representative households whose EHWT have been equipped with sensors (see Fig. 7.5). The characteristics of the EHWT (heat coefficient, power of the heating elements) are taken from representative products available on the market. The restrictive time intervals have been defined using the distribution of peak hours in France. The 18 distinct peak hours in France each correspond to a distinct proportion of the households. Finally, 8 realistic reference objective load curves have been considered.

Simulation results

Simulations have been conducted on the previously described dataset. The exact objective load curve have been employed but are omitted for confidentiality reasons. For illustration, examples for two fictional but representative different objective load curves are reported in Fig. 7.6. To quantify optimality losses, given the obtained load curve f b (t), we propose two following adimensional indexes, respectively corresponding to L 1 and L 2 norms.

q 1 = t f t 0 |f b (s) -f o (s)|ds t f t 0 |f o (s)|ds , q 2 = t f t 0 (f b (s) -f o (s)) 2 ds t f t 0 (f o (s)) 2 ds .
Results for the various reference objective load curves (with the best results among various probability law) and various values of n (number of tanks) and p (number of time-steps) are reported in Table 7.1 and Table 7.2. The presented heuristic has been implemented in Matlab R15a, and run on a Intel Core i7 (3.3 GHz) with 16 GB of RAM, using a single core and 25 MB of memory.

The simulations show satisfactory results and highlight the relevance of favoring diversity during the rescheduling procedure. Indeed, results with a high number of tanks have better optimality index that the one of the small instances. On the contrary, quality of the results is not increasing with the number of time-steps (when the number of tanks is too small)2 .

Interestingly, the proposed heuristic is computed in approx. linear time (with respect to the number of EHWT), and can be applied to very large sets of tanks. However, it is to be noticed that efficiency of a given probability law, and quality of the results depend on the shape of the objective curve and the flexibility yielded by the duration distribution. In our case, with realistic sets of tanks and reference objective curve, the flexibility is sufficient, except for the objective load curve 6 (see Table 7.2) in which steep increases and decreases of the consumption make optimization difficult. An example of the occurring phenomena (the actual curve being omitted for confidentiality reasons) is pictured in Fig. 7.7. This type of objective curve may require to split heating durations into two disjoint periods, to bring some additional flexibility. Some trial-and-error procedure may be necessary to address other problems, but this task has not been particularly tedious in the cases reported in the thesis.

Time (h)

Optimization with uncertainty

Formulation of the problem when only a fraction of EHWT are smart

The first problem (7.6) adresses the case when all tanks are smart so that they transmit information about their own state. In this section, we extend the problem to the case when only a fraction of the EHWT are smart3 . For each of the non-smart tanks, we assume that we can use an unbiased estimator êi 0 of e i (t 0 ) with standard deviation σ i e , and an unbiased estimator τ i 0 of τ i (t 0 ) with standard deviation σ i τ . We use the standard deviation to strengthen constraints (7.8) and (7.9) as

t 0 + ∆t i b ≥ t i 0 , t 0 + ∆t i b + d i b + 2σ i e ≤ t i f , (7.10) and ∆t i b + exp(k i ∆t i b )(τ i 0 + 2σ i τ ) ≤ t c . (7.11)
The rationale behind these constraints is that unfavorable cases when EHWT are heated during peak-hour or when not enough water has been heated at time t c are unlikely. Then, the problem considering uncertainty has the same formulation as (7.6), but the constraints are (7.7), (7.8) and (7.9) for smart and remotely controllable tanks, and (7.7), (7.10) and (7.11) for tanks that are only remotely controllable.

Simulations

Simulations have been conducted under the same conditions as in Section 7.3, with a 0% ratio of "smart EHWT", and with a Gaussian distribution centered about the estimator value of the tank energy, and with a (relatively large) standard deviation equal to one fourth of the estimator value. An example is shown in Fig. 7.8 and quantitative results are reported in Simulation results show a good level of robustness, using either criteria q 1 and q 2 . Analysis of Fig. 7.8 reveals that a spreading effect takes place in the right edge of the load curve. This phenomena, that looks like heat dynamics will be the subject of future works.

Summary

In this chapter, we have formulated an optimization problem and proposed a resolution method in the form of a heuristic for the optimal rescheduling of heating of large sets of EHWT. Numerical experiments conducted on real data stress the relevance of this method for parameters corresponding to French houses. The method is only marginally sub-optimal and appears to scale linearly with the dimension of the problem.

Similar studies should focus on parameters for other countries. Depending on the load curve, which relates to electricity producers constraints, the heuristic may require further developments. A main tuning parameters is the choice of the probability density employed in the iterative scheduling procedure. In particular, if sharp transients are to be considered in the load curve, the choice of the probability function may require further investigations.

On the numerical side, the approach could benefit from various classic techniques. In fact, parallelization, semi-lumping are possible ways to explore. This could help speed up the method, which is already reasonably fast, but can be important for large instances. A straightforward implementation of the presented methodology treats a representative set of one million EHWT in 6000 seconds. A satisfactory rescheduling is obtained after 3 random runs.

Chapter 8

Large-scale groups: modeling populations of EHWT with Fokker-Planck equations

Modélisation de grande populations de ballons en utilisant les équations de Fokker-Planck. Dans ce chapitre, nous proposons un modèle de comportement d'un parc (grand nombre) de ballons sous la forme d'un système d'équations aux dérivées partielles. Ces équations proviennent du comportement stochastique individuel de chaque ballon (décrit ici), qui permettent de décrire le comportement global de la population à l'aide des équations de Fokker-Planck, lorsque le nombre de ballons tend vers l'infini.

The approaches developed in Chapter 6 and 7 carefully coordinate the individual controls of the tanks. The coordinated decisions are to be sent to the automation devices embedded in every tank in the group. Leaving out laboratory experiments, this might seem unrealistic, since almost no existing EHWT is equipped with such automation devices. In the future, "smart EHWT"will be deployed and will gradually represent a non negligible part of the market but it will be a long time until every EHWT is smart. Interestingly, this is not really necessary. One can use only a subset of the whole group and expect the optimal strategy to simply smooth out the load curve, in a spirit similar to the observation made in Chapter 7 where it was stressed that for large number of tanks, the smoothing of the load curve induced by stochastic individual behavior is sufficient to solve the optimization problem at stake.

In this chapter, we model the behavior of the whole group of EHWT as a system submitted to external forcing terms. In the limit, we obtain a population dynamics. This approach can be used to test various kinds of strategies. For instance, effect of new schedule for time-of-use pricing can be quantified. Other perspectives are possible. They will be sketched in the conclusion of the thesis.

Interestingly, it is possible to recast the heating policies of a large set of EHWT into a distributed parameter framework. Following the works of [START_REF] Malhamé | Electric Load Model Synthesis by Diffusion Approximation of a High-Order Hybrid-State Stochastic System[END_REF] (recently used by [START_REF] Moura | Modeling Heterogeneous Populations of Thermostatically Controlled Loads Using Diffusion-Advection PDEs[END_REF]), we consider that the local (individual) control variables of EHWT are each defined according to stochastic processes. Then, we combine i) this randomness, ii) the diversity in the distribution of the states of the EHWT, iii) the randomness of the water consumptions, and we develop a PDE for a large group of EHWT. Classically, this takes the form of Fokker-Planck equations (see [START_REF] Risken | The Fokker-Planck Equation[END_REF]) governing the probability distributions of the population of EHWT. The work of Malhamé and Chong was originally focused on a mitigated load represented by a single state. For smart piloting applications under consideration here, extensions are necessary. The result are a rich system of PDE, which constitutes the main contribution of this chapter.

The chapter is organized as follows. In the group of EHWT, a single EHWT is a macroscopic but small compared to the whole group subsystem described by three state variables. Notations are given in Section 8.1. To account for the randomness of water consumption, we propose to represent single EHWT using a Markovian stochastic process in Section 8.2. Then, we introduce probability density functions of the population of EHWT and derive the Fokker-Planck equations in Section 8.3. A summary of the obtained input-output description of the EHWT group is reported in Section 8.4. Conclusions and perspectives are given in Section 8.5.

Model statement

Domains of definition of state variables

We recall (6.1), λ = Tcom-T in Tmax-T in . Let e max be the maximal energy that can be contained in the tank under the temperature T max . Then, by definition, the states a, τ , µ are subject to the following inequalities :

0 ≤ a ≤ e max , 0 ≤ τ ≤ λe max , 0 ≤ µ ≤ λe max , a + 1 λ (τ + µ) ≤ e max ,
Ω 0 ={(a, τ, µ)|a, τ, µ > 0, λe max < a + τ + µ, a + 1 λ (τ + µ) < e max } F 1 = Ω0 ∩ {(a, τ, µ)|µ = 0} F 2 = Ω0 ∩ {(a, τ, µ)|τ = 0} F 3 = Ω0 ∩ {(a, τ, µ)|λe max = a + τ + µ} F 4 = Ω0 ∩ {(a, τ, µ)|a + 1 λ (τ + µ) = e max }.
The following edges and vertices are considered

E 1 =F 1 ∩ F 2 E 2 =F 3 ∩ F 4 V 1 ={(λe max , 0, 0)} V 2 ={(e max , 0, 0)} V 3 ={(0, λe max , 0)} V 4 ={(0, 0, λe max )}.
In practice, z = (a, τ, µ) can only belong to E 2 (low energy edge, e.g. Fig. 5.2 (f)), Ω 0 (medium energy domain, e.g. Fig. 5.2 (a), (b) and (e)) and E 1 (high energy edge, e.g. Fig. 5.2 (c) and (d)). Faces F 1 to F 4 , and vertices V 1 to V 4 constitute boundaries of these three domains. Note that uniformly cold tanks can also stack in V 3 .

In the following, we note Ω = Ω 0 ∪ E 1 ∪ E 2 ∪ V 3 . Any index i = 0, 1, 2, 3 will refer to these sub-domains, respectively.

λe max V 1 a τ µ e max λe max λe max 0 F 3 Ω 0 F 2 E 1 F 1 F 4 E 2 V 3 V 4 V 2 Figure 8

Transient behavior

The transient behavior of z is driven by the physical phenomena in the tank. The dynamics is similar to the discrete one described in Section 6.1.2. Heating mostly induces a continuous variation of z. If τ, µ > 0 (i.e. z ∈ Ω 0 ), then the heating yields a decrease of τ and an increase of µ (see Fig. 5.2 (b)). Under certain conditions, a threshold effect can be observed: when τ reaches 0 (i.e. when z reaches F 2 ), then suddenly, all the unavailable energy µ becomes available, µ takes the value 0 and all its energy is transferred into a. This effect constitutes the evolution from Fig. 5.2 (b) to (c). This induces a discontinuity, transferring z to E 1 in which the heating has again a continuous effect on z, increasing the available energy a (see Fig. 5.2 (c) and (d)).

Heat losses also mostly induce a continuous variation of z, during which a and µ decrease while τ increases. The reverse threshold effect can also be observed: when z reaches F 3 , an entire layer of water goes under the temperature T com , which causes a to take the value 0 and z to jump in E 2 in which it will again vary continuously. This case can, for instance, be encountered with the rest of the profile displayed in Fig. 5.2 (d): when the uniform temperature in the tank reaches T com , a takes the value 0.

The drains have a very short duration compared to heating and heat losses. They yield quasi-instantaneous decrease of a and µ and increase of τ , having time of occurrence and magnitude depending on the user's habits. We model them as jumps only. Very large decreases can cause z to jump from E 1 to Ω 0 (or even to E 2 ), or from Ω 0 to E 2 .

A summary of jump dynamics is pictured in Fig. 8.2.

EHWT as a hybrid-state stochastic process

Stochastic process representation

The times of occurrence and magnitude of the drains are user-dependent. We model this dependency under the form of a stochastic process. This results in a model similar in spirit to the one-dimensional representation originally proposed by [START_REF] Malhamé | Electric Load Model Synthesis by Diffusion Approximation of a High-Order Hybrid-State Stochastic System[END_REF].

Heat losses and heating modeling as a drift

Under the assumption that the ambient temperature T a is equal to T in , the heat losses per unit of time for a t are equal tok Sρcp a t , where S is the cross-section of the tank, ρ and c p are the density and heat capacity of water, and k is the heat losses coefficient of the tank per unit of height (see again Fig. 5.2).

Likewise, for µ t , the heat losses is equal tok Sρcp µ t . The heat losses generate a positive effect on τ t which increases with the rate k Sρcp µ t . One has

v r (Z t ) = k Sρc p    -1 0 0 0 -1 0 0 1 0    Z t .
To this heat losses drift we add a drift v h (Z t ) due to power injection for tanks subjected to heating. In Ω 0 and E 2 , power injection p lowers τ t and increases µ t , so that one has v h (Z t ) = p 0 -1 1 T . In E 1 , the injected power only affects a t , and therefore v h (Z t ) = p 0 0 T . Finally, when the tank is heating,

v(Z t ) = v r (Z t ) + v h (Z t )
and otherwise v(Z t ) = v r (Z t ).

Drain as a Poisson process

The drains appear as a sequence of quasi-instantaneous events of various magnitudes. For this reason, we choose to model them as a non homogeneous compound Poisson process J t (see e.g. [START_REF] Applebaum | Lévy Processes and Stochastic Calculus[END_REF]). This assumption is questioned by the results in Chapter 4 that highlight an autocorrelation in the DHW consumption. However, the daily pattern represented by the mean consumption shows similarities with the intensity function of a Poisson process, and we assume that the temporal correlation fades out as the number of tanks grows. As a consequence, the time between jumps follows an exponential law of parameter χ(t), and the magnitude of jumps is characterized by a probability density function ω : Ω 2 × R + → R + . In words, a jump from the state z at time t maps to the rest of Ω with a distribution characterized by ω(z, •, t).

EHWT variability as a Wiener process

Some phenomena are not taken into account in the description above (e.g. diffusion). We choose to lump these into into a standard deviation parameter σ(Z t ) acting through a Wiener process. However, the closer Z t is to the boundary of Ω, the smaller the uncertainty should be on some of variables. Thus, certain components of σ(Z t ) must vanish at the boundaries of domains. Additionally, the heating increases uncertainty on the dynamics. In summary, the standard deviation takes the form σ(Z t ) = σ r (Z t ) + σ h (Z t ) for heating tanks, and σ(Z t ) = σ r (Z t ) otherwise.

On/off heating switch modeled with a Poisson process

The switchings between the two statuses h and r constitute a sequence of instantaneous events that can be controlled. We choose to also model it with a Poisson process of intensity α(Z t , t, S t ) = α St (Z t , t) (indexed on S t for transition from S t to the opposite one). This means that instead of exactly setting the switching times, two functions α r , α h define a probability to switch from one status to another, depending on the state Z t and time t.

Jump map: hybrid system modeling of the domain switch

The threshold jumps during heating (and in theory possibly during the rest phases) constitute hybrid deterministic jumps. When reaching a certain boundary, it maps a domain to another, depending on the status. The transition, in the framework of [START_REF] Goebel | Hybrid Dynamical Systems: Modeling, Stability, and Robustness[END_REF], gives z + (the value after the jump) as a function of z (the value before the jump). For the sake of clarity, a summary is given in Table 8.2. The maximal energy e max that can be contained in the tank is reached at point V 2 , when heating. We assume that the heating automatically switches off at this point for security reasons, which is characterized by a hybrid jump for the status from h to r. This jump is also presented in the jump map.

Fokker-Planck PDE for a large group of EHWT

EHWT group population representation

Representing a large group of tanks each having a 3-dimensional state leads to an unnecessarily large finite-dimensional system, which can be difficult to design controllers for. Rather, a probability density functions representation can be employed.

The main idea is to define 7 functions f r 0 , f h 0 , f r 1 , f h 1 , f r 2 , f h 2 , f r 3 (one for each domain Ω 0 , E 1 , E 2 and V 3 , one for each status r or h at the exception of V 3 in which only resting tanks can stack) which represent the population density of the tanks of a given status in a certain domain. These positive functions are subject to the balance

Ω 0 (f h 0 + f r 0 ) + E 1 (f h 1 + f r 1 ) + E 2 (f h 2 + f r 2 ) + f r 3 = 1.
The dynamics governing these probability functions are obtained from the preceding dynamics. We now detail them.

Fokker-Planck equation for a stochastic process

Population distribution can often be studied through the Fokker-Planck equation (see [START_REF] Risken | The Fokker-Planck Equation[END_REF]). For a set of independent Markov process in a state space Ω following the same generic stochastic equation where v(Z t , t), σ(Z t , t) ∈ R 6 , dW t is the integration with respect to a one-dimensional Wiener process W t , and dJ t is the integration with respect to a compound Poisson process of intensity χ(t) and whose compound distribution is represented with the probability density function ω (i.e. when a jump occurs on state z, the probability density function of transition to state y is represented with ω(z, y, t) at time t), the probability density function

dZ t = v(
In

Ω 0 E 1 E 2 If reaches F 2 × {h} F 3 × {r, h} V 1 × {r} V 2 × {h} V 4 × {h} Jump to E 1 × {h} E 2 × {r, h} E 2 × {r} V 2 × {r} E 1 × {h}
According to the transition

   a + τ + µ +    =    a + µ τ 0       a + τ + µ +    =    0 τ a + µ       a + τ + µ +    =    0 τ a + µ    s + = r    a + τ + µ +    =    a + µ τ 0   

Detailed expressions for the dynamics

Several observations can be made on the stochastic model presented in Section 8.2. First, for each tank, the stochastic process Z t defined by (8.2) constitutes a Markov process. Moreover, given a group of tanks, the independence of the stochastic process of each tank appears as a reasonable assumption, given that hot water consumptions of distinct households are usually not related. Therefore, we follow the work of [START_REF] Malhamé | Electric Load Model Synthesis by Diffusion Approximation of a High-Order Hybrid-State Stochastic System[END_REF], and derive the Fokker-Planck equations. In our case, for each domain and each status, this equation takes the form of a parabolic PDE. The hybrid nature of the stochastic process appears in the boundary conditions and yields an additional integral source term.

On each domain i = 0, 1, 2 (for Ω 0 , E 1 , and E 2 ), f r i and f h i are driven by a system of the form

∂ t f r i + ∇ z • [v r f r i ] = ∇ z • [D r i ∇ z f r i ] -(α r i + χ)f r i + α h i • f h i + S r i (f (•, t), z, t) ∂ t f h i + ∇ z • [(v r + v h i )f h i ] = ∇ z • [(D r i + D h i )∇ z f h i ] -(α h i + χ)f h i + α r i • f r i + S h i (f (•, t), z, t)       
(8.6) while f r 3 follows an ODE that will be stated later in (8.7).

Definition of parameters in each domain

Each term in the stochastic differential equation has a matching term in the partial differential equation (see e.g. [START_REF] Ethier | Markov Processes: Characterization and Convergence[END_REF][START_REF] Sato | Lévy Processes and Infinitely Divisible Distributions[END_REF]). For the sake of clarity, in Table 8.1, we give each stochastic term under consideration and its corresponding term in the PDE. The necessary steps of computations are omitted for brevity. The heat losses drift in the PDE has the same form as the one in the stochastic equation, i.e. In each domain, the drift v h i ∈ R 3 , the diffusion D r i (z, t), D h i (z, t) ∈ R 3×3 + , and the source terms S s i (f (•, t), z, t) ∈ R + for s = r, h have to be defined. The term caused by the Poisson process leads to source terms of various integral forms, depending on the probability density function ω i on each domain (e.g. ω 0 = ω |Ω 0 ). Source terms can also appear due to hybrid transfer from other domains in the form of an integral flow. For this reason, on F 2 and F 3 , we introduce the functions η 2 and η 3 defined as

E 1 p    1 0 0    d 1 (z)   
1 0 0 0 0 0 0 0 0

   S r 1 = 0 S h 1 = F 2 η 2 (y, z)(v 0 (y) + v h 0 )f h 0 (y, t)dy E 2 p    0 -1 1    1 2 d 2 (z)    0 0 0 0 1 -1 0 1 1    S r 2 = χ(

Domain Boundary Boundary condition

Ω 0 F 1 f r 0 (z, t) = 0 and f h 0 (z, t) = 0 F 2 f r 0 (z, t) = 0 and f h 0 (z, t) free F 3 f r 0 (z, t) free and f h 0 (z, t) free F 4 f r 0 (z, t) = 0 and f h 0 (z, t) = 0 They constitute valuable performance indexes for the system. In summary, the system admits the input-output representation depicted in Fig. 8.3. A natural goal is then to design controls (α(•, t) : Ω → R 6 + ) s.t. CB is as low as possible, while the total power demand P tot follows a given objective function.

E 1 V 1 f r 1 (z, t) free and f h 1 (z, t) = f h 2 (V 4 , t) V 2 f r 1 (z, t) = f h 1 (V 2 ,
As a mean of illustration, the probability density functions on E 1 , E 2 , and a representative segment of Ω 0 from the middle of E 2 to V 2 are shown in Fig. 8.4. Two profiles are shown. A fictional initial one, and the one subsequent to the following heating policy. We choose to promote heating (i.e. α r i high and α h i low) on E 1 , E 2 , and V 3 , and let the tanks rest (i.e. α r i low and α h i high) on Ω 0 . The profile varies as is shown in the figure, and tends to spread due to diffusion and integral drains. After some time, due to diffusion effects, a stationary profile representative of the cycle E 2 → E 1 → Ω 0 → E 2 should take place. Numerical treatment of the Fokker-Planck equations derived in this chapter will be the subject of further works.

Summary

In this chapter, we have explained the derivation of a model for a large group of EHWT. The input is a parameter defining the stochastic process of heating of each individual EHWT in the group. The outputs are the overall comfort variable defined in (8.8) and the total power demand (8.9). The dynamics are a collection of Fokker-Planck partial differential equations.

V 3 V 4 f 2 V 3 V 4 f 2 V 1 V 2 f 1 middle of E 2 V 2 f 0 V 1 V 2 f 1 middle of V 2 V 2 f 0 to E 1 from E 2 to Ω 0 to E 2 from V 3 to V 3

Part III

Conclusions and perspectives

Chapter 9

Conclusions and perspectives

Conclusions et perspectives. Dans cette partie, nous résumons les travaux effectués jusqu'ici en les mettant en relation. Les architectures de contrôle possibles pour chacune des méthodes abordées sont présentées. Enfin, des perspectives d'amélioration sont décrites.

In this thesis we have studied two central questions that we believe could be of interest for the development of prototypes of "smart EHWT". These questions focus on modeling and control.

On the modeling side, the model obtained in Chapter 3 has some interesting properties: it is computationally light and relatively accurate. Together with the domestic water consumption model of Chapter 4, it constitutes the basis for the input-output representation of Chapter 5, which could be implemented in a "smart EHWT".

On the optimal control side, we have considered three typical use-cases, sorted by ascending sizes. Very generally, several possible decisional architecture can be considered; they are pictured in Fig. 9.1. For small size the proposed reformulation as a MILP/MIQP yields to relatively light computational efforts, so that embedded resolution of optimal strategies is doable. This case requires that all tanks should be smart (see Fig. 9.2), and corresponds to cases (b) and (c) in Fig. 9.1. For medium size, it has been shown that even with large dimension, the problem of optimal scheduling can be solved with little sub-optimality thanks to a carefully designed heuristic. This case only requires that the tanks are controllable (see cases (a) and (b) in Fig. 9.1). When the dimension gets sufficiently large, the limit case can be approached by a partial differential equation, of the Fokker-Planck type, which opens new perspectives.

We now wish to draw some perspectives on the presented work.

1. The modeling efforts could, at little expense, be generalized to the case of thermodynamics and solar tanks, which are important systems in the context of DSM. However, it is likely that no plateau will be observed, despite the positive impact of stratification. The culprit is the spatially distributed nature of the internal heat exchangers. For this reason, no straightforward derivation of a multi-period model can be determined at this stage. However, some similar simplifications could be done, but investigations are needed.

2. The MILP/MIQP approach developed in this thesis yields good results but is strongly limited by the dimensionality of the problem. At this stage, we cannot address groups

Figure 1 . 1 :Figure 1

 111 Figure 1.1: Simplified scheme of an EHWT. The temperature sensors set for experimentation have been depicted. These are not found on off-the-shelf EHWT.

Figure 1 . 3 :

 13 Figure 1.3: Thermal image of an EHWT (obtained with an infrared camera).

Figure 1

 1 Figure 1.4: Input-output representation of the tank.

Figure 1 . 5 :Figure 1

 151 Figure 1.5: Schéma simplifié d'un CEJ. Les capteurs de température ajoutés pour les mesures expérimentale ont été représentés, ils ne sont pas présents sur les CEJ du commerce.

Figure 1 . 7 :

 17 Figure 1.7: Image thermique d'un CEJ (obtenue à l'aide d'une caméra infrarouge).

Figure 1 . 8 :

 18 Figure 1.8: Représentation entrée-sortie d'un CEJ.

Figure 2

 2 Figure 2.1: Plumes of turbulent natural convection over an exchanger, from [Bla10].

Figure 2 . 2 :

 22 Figure 2.2: Variations of the temperature profile during a heating period (a) and during a draining period (b) (blue: model prediction, black: experimental values).

Figure 3

 3 Figure 3.1: Definition of the timeline.

Figure 3

 3 Figure 3.2: Energy balance for the integral form heating model.

Figure 3

 3 Figure 3.3: The Stefan problem (with moving boundary) for modeling draining.

  Figure 3.4: Comparisons of the physics-based model in § 2.1.4 (a) and the proposed multi-period model (b), against experimental data (black), simulation are in blue.

Figure 4 . 1 :

 41 Figure 4.1: An example of sequence of drains for an household over 48 h. (experimental results)

Figure 4 . 2 :

 42 Figure 4.2: Experimental probability distribution of drain magnitudes (over the whole recorded sequences). The smoothed plot highlights 4 distinct peaks. (experimental results)

Figure 4

 4 Figure 4.3: Two examples of averaged daily cumulative number of drains. Experimental results are in blue, and and smoothed version is in orange.

Figure 4 Figure 4

 44 Figure 4.5: Experimental CDF for the adjusted increments compared to the Weibull and exponential CDF.

Figure 4 . 7 :Figure 4 Figure 4 . 9 :

 47449 Figure 4.7: Probability plot of the WACD simulation increments against experimental increments, for the whole set of data. (our proposed model)

Figure 5 . 1 :

 51 Figure 5.1: Temperature profile inside an EHWT after a drain. Available, delay and reserve energies correspond to area defined by the profile.

T

  -T (y))dy µ = Sρc p yc 0 (y)dy where ρ and c p are the density and the heat capacity of water, respectively, and y c is defined as y c min{y|T (y) = T com }.

Figure 5 . 2 :

 52 Figure 5.2: Schematic variation of the state z during heating and draining.

Figure 5

 5 Figure 5.3: Available, delay and reserve energies and schematic reconstruction of the original temperature profile.

Figure 5 . 4 :Figure 5 . 5 :

 5455 Figure5.4: Variation of z = (a, τ, µ) over a 24 h run. a: available energy, τ : delay energy, µ: reserve energy. Blue: experimental data. Red: simulation using the multi-period model. Black: Simulation using a single-zone model (with reconstructed variables).

Figure 6

 6 Figure 6.1: Strengthening of (I t ).

  red), in % of max power, and normalized prices c (black) green), τ (blue), µ (black), in % of tank capacity

Figure 6 . 2 :

 62 Figure 6.2: Optimization for a two price scenario.

  red), in % of max power, and normalized prices c (black) green), τ (blue), µ (black), in % of tank capacity

Figure 6 . 3 :

 63 Figure 6.3: Optimization for a spot prices scenario.

  in % of tank capacity

Figure 6 . 4 :

 64 Figure 6.4: Simultaneous optimization of four tanks.

Figure 7

 7 Figure 7.1: Day-ahead SPOT market price of electricity in France on 2/9/16. The shaded regions corresponds to the peak consumption of EHWT in France. The zone of lowest prices lies out of this region.

Figure 7 . 2 :

 72 Figure 7.2: Heating starting times and durations.

Figure 7

 7 Figure 7.3: Detailed process of the heuristic. (a) at the start of an iteration, a cumulative curve has been determined by scheduling some of the heatings. (b) a new heating is to be scheduled, considering its duration, some intervals of possible starting times are determined from the residual curve. (c) along these intervals, probability density functions are defined, they favor the boundaries of the intervals. (d) randomly, a starting time is selected and a new cumulative curve is computed.

Figure 7 . 4 :

 74 Figure 7.4: Example of three possible parabolic distribution used in practice on the admissible intervals.

Figure 7 . 5 :

 75 Figure 7.5: Typical distribution of durations derived from a representative panel of users, (normalized scale).

Figure 7

 7 Figure 7.6: Two examples of objective load curve (blue) and final realisations (orange). (a) n=500, p=100 (b) n=5000, p=1000.

Figure 7 . 7 :

 77 Figure 7.7: Unsuccessful scheduling for steep increase and decrease of the consumption. n=5000, p=1000.

  λe max ≤ a + τ + µ from which we define Ω 0 the following open polyhedron of R 3 and its faces (see Fig. 8.1):

Domain Drift v hi

  Diffusion D r i (z) and D h i (z) Source terms S r i (f (•), z, t) and S h i (f (•), z, t)

Figure 8 . 4 :

 84 Figure 8.4: Variation of probability density function, given a heating policy (red: heating, blue: resting).

The context of demand side management

  T pt , T px , T p∆ functions for the temperature of the plateau

							Contents Contents
	notation Acronyms notation meaning meaning		unit	unit
		a A		available energy forbidden line		-	J
		c p c		specific heat capacity cost function		W -1	J•kg -2 K -1
		C s C i		space-time domain for the Stefan equation total energy consumption of tank i on the horizon J	-
	DHW C ACD autoregressive conditional duration total domestic hot water consumption CB population braking the comfort constraints		-	m 3
	F I , F II , F III multi-period mappings for each phase F i filtration at the drain i CDF cumulative distribution function CEJ chauffe-eau Joule d i heating duration of tank i d t total energy drain at time period t Chapter 1		s J	--
	g DHW domestic hot water gravitational acceleration D diffusion		-	m•s -2
	h I k DSM demand side management height of the tank time interval heat losses to the ambient coefficient DWC domestic water consumption EACD Exponential ACD e i energy in tank i êi 0 estimator of the initial energy in tank i e p minimum energy security margin Introduction		J J J	m -s -1
	l EHWT electric hot water tank function for the plateau height dynamics e max maximum energy in the tank		J	m•s -1 W -1
	L c LP e i max	characteristic length in the tank linear program maximum energy in the tank i		J	m
	M M j MILP mixed integer linear program mean number of drains magnitude of the drain number j MIQP mixed integer quadratic program E edge f o objective load curve 1.1		-W	-m 3
	N k ODE f a	number of drain for day k ordinary differential equation initial load curve		W	-
	P W PDE f b	rescaled power injection per unit of length partial differential equation final load curve		W	K•s -1
	Ra QP f o	Rayleygh number quadratic programs objective load curve		W	-
	Ri SOS f i r	Richardson number special ordered set residual load curve at step i		W	-
		s f r i		function for the plateau height dynamics resting population density function		-	m•s -1
		S f h i		cross-section of the tank heating population density function		-	m 2
	t 0 , t f F Greek letters T p	initial and time temperature in the tank tuning parameter face		--	s o C
		∆T P tot		temperature spread in the tank total power demand		W	o C
		T 0 P 1		initial temperature profile cost minimization problem		-	o C
		T a P 2		temperature of the ambient objective load curve problem		-	o C
	notation α d T max T com P Notations relative to the EHWT groups meaning turbulent diffusion maximum temperature in the tank (set by the user) o C unit comfort temperature (set by the user) o C probability law on S i -q 1 , q 2 optimality index for L 1 and L 2 norms -m 2 •s -1 α th thermal diffusion T in inlet water temperature o C S source terms in the Fokker-Planck representation -m 2 •s -1 α V volumetric coefficient of thermal expansion of water T s temperature at the surface of the heating element o C S i set of admissible starting times -K -1 T ∞ characteristic temperature of the water in the tank o C S t tank status at time t (heating or at rest) -
	δ d Greek letters parameter of a Weibull distribution turbulent to thermal diffusion ratio T p temperature of the plateau t i 0 , t i beginning and ending horizon time for tank i f i standard distribution for the increments i consumption time t c ζ tuning parameter µ reserve energy ν kinematic viscosity of water T b temperature of the Stefan homogeneous zone T min tuning parameter T u , T v , T discontinuity point sequence u heat injection in the tank notation meaning u i power control function for tank i unit α s i maximal power injection u max coupling control term from status s, in domain i -α tuning parameter u i heating power of tank i max -U admissible power control set	s s W W W -	----J o C m 2 •s -1 o C o C -W
	ξ θ β ∆t	U h v v t , w t V	tuning parameter Heavyside function overall heat transfer coefficient tuning parameter tuning parameter heating starting time power injection in a and µ at time period t vertex	-s	J -	--W•m -2 K -1 m 1-ξ •s -1 K -ζξ
	κ η	v d v	scale parameter of a Weibull distribution drain velocity jump distribution at boundaries drift term	-	-	-m•s -1
	Π ρ λ τ υ φ i φ Φ χ φ t σ σ i e σ i τ 0 τ i	v nc Y y c y v max x x p x b Y z Z	tank perimeter water density natural convection velocity comfort to maximum energy in the tank ratio internal state for MILP/MIQP delay energy normalized increments between two successive drains -mean conditional duration for the normalized increment i -m m•s -1 -kg•m -3 J -tuning parameter exchange coefficient minimal height at comfort temperature m intensity of the Poisson process -s -1 time interval between two successive drains s energy flow from µ to a at time period t J m -1 tuning parameter m•s -1 height in the tank m height of the plateau m height of the Stefan homogeneous zone m standard deviation term in stochastic dynamics -admissible internal state for MILP/MIQP -standard deviation on the estimator êi 0 J three-dimensional energy state -standard deviation on the estimator τ i 0 J estimator of the initial delay energy in tank i J stochastic state process -
	ω 0 , ω 1 , ω 2 ACD parameters ω jump distribution	-		-
	Ω Ω		space domain state domain	-		-

o C

Part I The EHWT : behaviour and proposed representation Chapter 2 A physics-based representation of EHWT Un modèle physique pour les CEJ. Dans

  ce chapitre, nous présentons un modèle de CEJ prenant la forme d'un système couplé de deux equations aux dérivées partielles, suivi de simulations et d'une validation expérimentale. Le modèle permet de reproduire le profil de température dans le ballon au cours du temps.

	[BMP15a] N. Beeker, P. Malisani, and N. Petit. A distributed parameters model for
	electric hot water tanks. In Proceedings of the American Control Conference,
	ACC, 2015.
	[BMP15b] N. Beeker, P. Malisani, and N. Petit. Dynamical modeling for electric hot
	water tanks. In Proceedings of the Conference on Modelling, Identification and
	Control of Nonlinear Systems, MICNON, 2015.
	[BMP16a] N. Beeker, P. Malisani, and N. Petit. Discrete-time optimal control of electric
	hot water tank. In Proceedings of the 11th IFAC Symposium on Dynamics and
	Process Systems, including Biosystems, DYCOPS, 2016.
	[BMP16b] N. Beeker, P. Malisani, and N. Petit. Modeling populations of electric hot
	water tanks with Fokker-Planck equations. In Proceedings of the 2nd IFAC
	Workshop on Control of Systems Governed by Partial Differential Equations,
	CPDE, 2016.
	[BMP16c] N. Beeker, P. Malisani, and N. Petit. An optimization algorithm for load-
	shifting of large sets of electric hot water tanks. In Proceedings of the 25th
	International Conference on Efficiency, Cost, Optimization, Simulation and
	Environmental Impact of Energy Systems, ECOS, 2016.
	[BMP16d] N. Beeker, P. Malisani, and N. Petit. Statistical properties of domestic hot
	water consumption. In Proceedings of the 12th REHVA World Congress CLIMA,
	2016.
	Liste des brevets
	[BMC15a] N. Beeker, P. Malisani, and A.S. Coince. Ballon Intelligent -1554897 (patent
	applied for), 2015.
	[BMC15b] N. Beeker, P. Malisani, and A.S. Coince. Ballon Intelligent -EDP -1554898
	(patent applied for), 2015.

Table 2 .

 2 three nearby elements permitting a power injection up to 2200 W. The dimensions of the water tank are specified in Table2.1. 1: Specifications of the EHWT used in experiments.

	Volume	L	200
	Length	m	1.37
	Maximal power	W	2200
	Heat losses coefficient W•m -2 K -1 0.66

To validate the model, experiments have been conducted in the facilities of EDF Lab Research Center, on an Atlantis ATLANTIC VMRSEL 200 L water tank. The power is injected via

  8 o C+ Time Distr. 53.9 % 22.9 % 10.7 % 5.1 % 7.4 % 2435.6 s Table 2.3: Comparison of absolute difference between experimental results and model predictions. Percentage of sample for each error interval.

  Table 3.1: Comparison of absolute difference between experimental value and model prediction; computational time. Percentage of sample for each error interval.

	8 o C+	Time
	Physic-based model 53.9 % 22.9 % 10.7 % 5.1 % 7.4 % 2435.6 s
	Multi-period model 82.3 % 14.5 % 2.0 % 0.6 % 0.6 %	4.6 s

  This new relation has the linear form required in linear programming formulations, and if e max is sufficiently large, only a small feasible regime (see Fig.6.1) is left out of the optimization problem.Given this new relation, we define the polytope

6.1). Instead of considering (I t ) we (conservatively) consider (I t ) which has the linear form M τ t ≥ µ t , (I t ) chosing M > 0.

Table 7

 7

	Number of tanks n Number of time-steps p	q 1	q 2	Computation time (s)
	500	500		0.0166 0.0162	1.8
	500	1000		0.0126 0.0122	2.1
	500	2000		0.0186 0.0179	4.0
	5000	500		0.0032 0.0031	17.6
	5000	1000		0.0028 0.0029	22.3
	5000	2000		0.0030 0.0032	33.9
	50000	500		0.0024 0.0021	131.5
	50000	1000		0.0018 0.0017	201.0
	Table 7.1: Results for objective load curve 1.
	Objective load curve	q 1	q 2	Computation time (s)
	1	0.0028 0.0029	22.3
	2	0.0036 0.0041	18.6
	3	0.0039 0.0042	18.5
	4	0.0037 0.0044	17.9
	5	0.0032 0.0030	25.6
	6	0.0233 0.0245	20.3
	7	0.0038 0.0042	15.4
	8	0.0049 0.0065	17.2

.2: Results for objective load curves 1 to 8 (5000 tanks, 1000 time-steps).

  Table 7.3.

	Objective load curve	q 1	q 2	Computation time (s)
	1	0.0338 0.0365	23.2
	2	0.0426 0.0446	17.8
	3	0.0336 0.0370	18.3
	4	0.0372 0.0376	16.9
	5	0.0319 0.0397	17.6
	6	0.0410 0.0609	20.6
	7	0.0335 0.0437	18.0
	8	0.0470 0.0544	16.9
	Table 7.3: Results for objective load curves 1 to 8 (5000 tanks, 1000 time-steps), with
	uncertainty.

  .1: Domains of definition of the state variables: F 1 to F 4 are the four faces of the open domain Ω 0 . V 1 to V 4 are its vertex, and E 1 and E 2 are two interesting edges. The tank state z belongs to Ω 0 , E 1 or E 2 .

			Drains (Poisson process)
			Heating (hybrid jumps)
			Heat loss (hybrid jumps)
	E 2	Ω 0	E 1
	Low energy		High energy
		Comfort limit

Figure 8.2: Transient between sub-domains.

  Z t , t)dt + σ(Z t , t)dW t + dJ t (8.3)

Table 8 .

 8 2: Jump map.(when the number of stochastic processes tends to infinity) is given by the Fokker-Planck equation (see e.g.[START_REF] Applebaum | Lévy Processes and Stochastic Calculus[END_REF])∂ t f (z, t) = -∇ z • [v(z, t)f (z, t)] + ∇ z • [D(z, t)∇ z f (z, t)]

	+ χ(t)		(f (y, t) -f (z, t))ω(z, y, t)dy	(8.4)
		Ω
	for (z, t) ∈ Ω × R + , where
	D(z, t) =	1 2	σ(z, t) • σ t (z, t) ∈ R 3×3 .	(8.5)

  t) Table 8.3: Definition of distributed parameter equation term on each domain.

	Ω 0 f r 0 f h Ω 0 (y)ω 0 (y)ω 0 (y, z, t)dy + χ(t) E 2 f r 2 (y)ω 2 (y, z, t)dy 0 (y, z, t)dy + F 3 η 3 (y, z)v 0 (y, t)dy 0 (y)f h 2 f h +χ(t) E 2 (y)ω 2 (y, z, t)dy
	S h 2 = χ(t)

  Figure 8.3: Input-output representation of the group of EHWT. and the total power demand is P tot (t) =

					t) and f h 1 (z, t) free
	E 2	V 3 V 4	f r 2 (z, t) free and (v h 2 + v r (z))f h 2 (z), t) = α r 2 (z), t)f r 3 (t) f r 2 (z, t) = 0 and f h 2 (z, t) free
		Table 8.4: Boundary conditions.
			PDE system	CB(t)
		α(•, t)	on f		P tot (t)
			2
			i=0	Ω	f h i (z, t)dz.	(8.9)

such as France's nuclear installation, which represented 76,3% of the country's electricity production in

[Bil15] 2 Similar effects have been observed, worldwide 3 the market share of electric heater is 35% in Canada, see[START_REF] Aguilar | Domestic water heat and water heater energy consumption in Canada[END_REF], 38% in the U.S, see[START_REF] Ryan | Water heater market profile[END_REF], 45% in France, see[START_REF] Msi | Marché des équipements domestiques de production d[END_REF] 4 Other examples are the Economy 7 and Economy 10 differential tariff provided by United Kingdom electricity

comme le parc nucléaire français, qui représentait 76,3% de la production nationale d'électricité en 2015[Bil15]

la part de marché du CEJ dans les moyens de chauffe est de 35% au Canada [AWR05], de 38% aux U.S.A [RLLL10], de 45% en France[START_REF] Msi | Marché des équipements domestiques de production d[END_REF]

les tarifications différentielles Economy 7 et Economy

au Royaume-Uni sont d'autres exemples typiques

[BMC15c] N. Beeker, P. Malisani, and A.S. Coince. Ballon Intelligent -Trois zones -

This model enrichment can be related to the numerics oriented works of[START_REF] Vrettos | Load frequency control by aggregations of thermally stratified electric water heaters[END_REF] who proposes to consider an additional natural convection term in the finite difference scheme discretizing a one-dimensional convection-diffusion equation.

The Richardson number (2.1) is a dimensionless number representing the relative importance of natural convection compared to forced convection[START_REF] Zurigat | Influence of inlet geometry on mixing in thermocline thermal energy storage[END_REF]

This correlation is influenced by the geometry of the inlet nozzle[START_REF] Ibrahim | Reviex of water-heating systems: General selection approach based on energy and environmental aspects[END_REF][START_REF] Zurigat | Influence of inlet geometry on mixing in thermocline thermal energy storage[END_REF]

usually in scheduling problems, the assignment of a job to a machine is required

In fact, a discrete time-mesh refinement brings significant additional performance only with an increasing number of tank. This point should require additional investigations.

By definition, all the tanks under consideration are remotely controllable

hybrid systems are dynamical systems that are subjected to discontinuities in their state in certain parts of their domain. A complete framework for hybrid system is presented in[START_REF] Goebel | Hybrid Dynamical Systems: Modeling, Stability, and Robustness[END_REF].

Remerciements

Time (h)

.8: Uncertain case: example of objective load curve (blue) and scheduling prevision (orange) and final realization when the duration of each tank is known (yellow). n=5000, p=1000.

In our stochastic model, each EHWT is defined by the state vector

and its heating status S t ∈ {r, h}, which is also of stochastic nature. From now-on, the indexes {r, h} refer to "rest" and "heating", respectively. The dynamics of Z t is governed by the phenomena described in Section 8.1.2. In each domain Ω 0 , E 1 , E 2 , the state vector Z t changes continuously (due to heating and heat losses) or discontinuously (due to drains), with respect to a stochastic differential equation, constituting the flow map, that will be stated below. A jump of hybrid nature 1 appears when Z t reaches F 2 , F 3 , V 1 or V 4 . This discontinuity, in accordance with the terminology of [START_REF] Goebel | Hybrid Dynamical Systems: Modeling, Stability, and Robustness[END_REF], constitutes the jump map. A correspondence between the phenomena and their stochastic/hybrid counterparts is detailed in Table 8.1. The stochastic differential equations of Z t and S t are

Flow map: Stochastic process dynamics on each domain

where, using the arbitrary choice we propose,

• v(Z t )dt is the drift component which represents the heat losses and/or heating effects;

• dJ t is the infinitesimal integration with respect to a 3-dimensional compound Poisson process J t representing the jump effects of drains on a t , τ t and µ t ;

• uncertainties are lumped into a standard deviation term σ(Z t ), integrated with respect to a Wiener process (or standard Brownian motion) W t ;

• dN t is the infinitesimal integration with respect to a one-dimensional Poisson process N t representing the status switch between h and r.

Expressions of v, σ and of the random characteristics of J t depend on the domain and the status of the EHWT. We now detail them and discuss the choices.

Chapter 8. Large-scale groups follows: η 2 (y, z) = 1 if a z τ z µ z T = a y + µ y τ y 0 T and 0 otherwise; and η 3 (y, z) = 1 if a z τ z µ z T = 0 τ y a y + µ y T and 0 otherwise.

Details are reported in Table 8.3. Finally, the exchange terms α s i (z, t) can be chosen as they are control-dependent, while χ(t) does not depend on space.

Boundary conditions

The domain Ω 0 has 4 boundaries (F 1 to F 4 ), while E 1 and E 2 have 2 boundaries, each in the form of vertices. Boundary conditions stem from the behavior of the stochastic process. A special case is the boundary V 3 where uniformly cold tanks stack.

On the boundaries of the domain, we have :

This allows to define boundary conditions of the Dirichlet or free boundary types, except in V 3 . Their definitions stem from exchange between the domains: free boundary corresponds to the case where tanks flowing outside the domain flow inside another domain as a source term (hybrid jumps). On the contrary, the zero Dirichlet boundary conditions correspond to the fact that no new tank can enter the system (the population is fixed). A summary is presented in Table 8.4. The vertex V 3 constitutes a free boundary for f r 2 , in which the population of completely cold tanks can stack. The population f r 3 (t) at this point can be heated and constitutes the input flow for f h 2 through the boundary condition

Remark 6. In (8.7), for ease of notation, the integration along E 2 is denoted with a 3-dimensional element dy to avoid the introduction of a parametrization function of E 2 . When no ambiguity is possible, similar notational simplifications are used in the rest of the chapter.

Input-output model

The input of the system is the set of functions (defined over space and time)

T which determine the intra-domain migration between the populations of heating and resting tanks.

To control the group of tanks, several indicators can be interesting. They are the output of the proposed model. Among them, the mass (number) of tanks breaking the comfort constraints is

Chapter 9. Conclusions and perspectives controllable EHWT decision center of EHWT having more than 4 elements. Fortunately, the constraints are decoupled and the resolution of single-EHWT problems is easy and fast. This naturally stresses that decomposition-coordination methods [START_REF] Cohen | Optimization by decomposition and coordination: A unified approach[END_REF] are well suited for this extension.

A promising type of decomposition is the "price-decomposition".

3. The heuristic proposed in this thesis can be the subject of numerous practical refinements. In particular, if the profiles to be attained are less smooth than the one presented in the thesis (e.g. featuring steep descents followed by sharp ascents), or if it is desired to address the problem on 48 h with two successive load curves, then more advanced heuristics are needed. In this direction, we have already developed refined strategies for cases of practical industrial interest. These remain out-of-the scope of the thesis. The questions at stake are mostly of computational nature.

4. In the case when not all the EHWT are "smart EHWT", the application of the heuristic for medium-scale groups lead to a subsidence (sinking) of the realized load curve compared to the desired one. The subsidence features a diffusion comparable to the heat equation. To attain the objective load curve, it may be useful to aim for a pre-processed load curve (e.g. with some sort of anti-diffusion) to sharpen the realized load curve. To address this problem, the effect of uncertainty on the load curve has to be quantified, formally.

5. At last, we believe that the Fokker-Planck approach has a lot of potential, even for real applications. The developments presented here are just a first step. We expect that exploiting the whole class of methods of the Fokker-Planck type can be a most effective way of representing and solving optimal control problems for the population of EHWT, partially blended with "smart EHWT", already installed in France. Our choice of modeling is still arbitrary, and numerous other possibilities must be explored. The next steps should address the control problems based on this input-output description. A question to be solved can be formulated as follows: how to design α(•, t) so that the power demand approaches some desirable profile while limiting or minimizing the discomfort? This problem belongs to the class of optimal control (tracking) of distributed parameter systems with in-domain actuation.

Résumé

Cette thèse s'intéresse au développement de stratégies de décalage de charge pouvant être appliquées à un parc de chauffe-eau Joule (CEJ).

On

Abstract

This thesis focuses on the development of advanced strategies for load shifting of large groups of electric hot water tanks (EHWT).

The first part of this thesis is dedicated to representing an EHWT as an input-output system. The idea is to design a simple, tractable and relatively accurate model that can be implemented inside a computing unit embedded in a "smart EHWT", for practical applications of optimization strategies. It includes in particular a phenomenological multi-period model of the temperature profile in the tank and a model for domestic hot water consumption.

The second part focuses on the design of control strategies for a group of tanks. Three use-cases are studied. The first one deals with a small number of smart and controllable EHWT for which we propose a discrete-time optimal resolution method. The second usecase addresses a medium-scale group of controllable tanks for which we propose a heuristic to optimally schedule the heating periods. Finally, we present the modeling of the behavior of an infinite population of tanks under the form of a Fokker-Planck equation.

Keywords

Electric hot water tank; Energy storage; Supply of hot water; Domestic water consumption; Dynamic Optimization; Multiperiod model; MILP; Heuristics; Fokker-Planck