Tadesse Getnet

Ayele

Keywords: District heating networks, Electricity distribution networks, Extended energy hub approach, Integrated load flow study, Integrated optimal power flow, Multi-carrier energy networks, Multi-energy system, Optimal placement, Thermal storage

Recent literature shows that there is a significant potential of decarbonisation and efficiency improvement that can be achieved through the synergy from multi-energy systems (MESs). Coupling technologies, such as co-generation plants, heat pumps and thermal storages are widely recommended as means of unlocking additional flexibility and increasing the penetration of renewables in the heating and electricity sectors. In view of that, the size and number of coupling technologies, such as combined heat and power plants and heat pumps (HPs), being installed in the heat distribution networks are increasing. As these technologies are exclusively managed by the district heating network operators, their operation sometimes becomes suboptimal from the electricity network point of view, and they (in particular large HPs) may cause overloading of the low voltage electricity distribution networks.

Integrated simulation and optimisation models are required to exploit the synergies effectively without compromising the constituent distribution networks of MES. Such models are not yet well developed. The conventional single-energy-carrier simulation tools are not capable of capturing key operating parameters of the multi-carrier distribution networks either.

A novel methodology for simulation and optimisation of MES is developed in this thesis based on an Extended Energy Hub (EEH) approach. The general framework is first developed in modular form so that it can be easily adapted for any type of multi-carrier energy networks. The framework is then used to develop the details of an integrated load flow model governing coupled heating and electricity distribution networks. Various load flow case studies with radial and meshed topologies are considered for demonstration and numerical validation of the proposed model.

The load flow model is further combined with a particle swarm optimisation algorithm in order to conduct integrated optimal power flow studies. Its contribution to the state of art is demonstrated by studying the optimal placement of coupling technologies, such as HPs and boilers in coupled heating and electricity distribution networks. The capacity of the model is further illustrated by exploiting the synergies using HPs together with thermal storage in the presence of intermittent renewables and variable electricity price signal.

It is shown that the EEH-based simulation and optimisation methodologies proposed in this thesis are very effective, flexible and easily scalable in capturing the key operating parameters of integrated electricity and district heating networks. The models can be used as a platform for further studies on integration of smart grids and smart thermal networks.

Résumé substantiel Préface

Cette thèse est basée sur une recherche réalisée dans le cadre du programme conjoint de doctorat Erasmus Mundus SELECT + 'Environomical Pathways for Sustainable Energy Services' coordonné par KTH Royal Institute of Technology, Stockholm, Suède. La recherche faisait partie du projet « Analyse des opportunités du réseau électrique décentralisé et de l'interaction du réseau de chauffage urbain ». Il a été financé avec le soutien de l'Agence exécutive pour l'éducation, l'audiovisuel et la culture (EACEA) (FPA-2012-0034) de la Commission Européenne.

La recherche se concentre sur l'intégration optimale des réseaux d'électricité et de chauffage urbain. Elle a été réalisée au sein de deux instituts (d'origine et d'accueil). Le département des systèmes énergétiques et de l'environnement (DSEE) de IMT Atlantique, Nantes, France (qui fait partie de l'UMR CNRS 6144 GEPEA) a été l'institut d'origine où un peu plus de la moitié de la période de recherche a été réalisée. Le reste du temps s'est déroulé au sein de la division Heat and Power Division (EKV) du Department of Energy Technology (EGI), à la School of Industrial Engineering and Management (ITM) de KTH Royal Institute of Technology, Stockholm, Suède, qui a servi d'hôte institut.

En résumé des recherches menées, cette thèse présente de nouvelles approches de modélisation, simulation et optimisation qui peuvent être utilisées pour exploiter les synergies du couplage des réseaux de chauffage urbains et d'électricité sans compromettre les réseaux de distribution respectifs.

« Exploiter les synergies des réseaux couplés de distribution d'électricité et de chaleur : Modélisation, simulation et optimisation basées sur une approche étendue de hub énergétique »

Contexte

L'énergie est un apport indispensable pour gérer les activités humaines quotidiennes. Tous les secteurs, de l'agriculture à l'industrie, de l'éducation au système de santé, du résidentiel au commercial, ont besoin d'une ou plusieurs formes d'énergie pour fonctionner correctement. La quantité d'énergie consommée au niveau des ménages est également l'un des indicateurs importants pour estimer la qualité de vie d'un pays [START_REF] Joyeux | Household energy consumption versus income and relative standard of living: A panel approach[END_REF][START_REF] Al-Mulali | Exploring the bi-directional long run relationship between energy consumption and life quality[END_REF][START_REF] Nadimi | Modeling of quality of life in terms of energy and electricity consumption[END_REF]. La consommation mondiale d'énergie continue de croître avec l'urbanisation et l'industrialisation. Il est rapporté que la consommation d'énergie dans le monde a augmenté de 2,3% en 2018, un taux qui est près du double du taux moyen depuis 2010 [START_REF]IEA -International Energy Agency[END_REF]. La croissance de l'économie mondiale n'est pas le seul facteur explicatifs de ce taux de croissance record. L'augmentation de la demande de chauffage et de refroidissement due au changement climatique, qui se manifeste par des températures moyennes d'hiver et d'été dépassant leurs records historiques dans certaines régions, a également contribué à cette augmentation. Bien que la capacité installée d'énergies renouvelables ait doublé pour la même période, le taux d'EnR n'étant pas suffisant pour faire face au taux de croissance de la demande énergétique, une partie importante de l'augmentation de la demande a donc été couverte par l'utilisation accrue de combustibles fossiles. En conséquence, les émissions de CO2 liées à l'énergie ont augmenté de 1,7%. Cette situation pourrait conduire à une boucle circulaire de cause à effet : plus d'émissions de CO2conditions météorologiques pires -croissance de la demande d'énergie -plus d'émissions de CO2.

Il est également signalé qu'en raison de l'augmentation de la demande mondiale d'énergie, d'une part, et de la rareté et des préoccupations environnementales des combustibles fossiles, d'autre part, l'énergie est devenue l'un des principaux défis auxquels l'humanité est confrontée au cours du siècle actuel [START_REF] Mohammadi | Energy hub: From a model to a concept -A review[END_REF]. Une transition énergétique vers un système énergétique plus renouvelable et durable est largement préconisée comme solution. Cette transition énergétique est une voie par laquelle le secteur mondial de l'énergie doit passer de « fossile » à « zéro carbone » d'ici la seconde moitié de ce siècle. On estime que jusqu'à 90% des objectifs de décarbonation peuvent être potentiellement couverts par l'application de mesures drastiques sur les énergies renouvelables et l'efficacité énergétique [START_REF]Energy Transition[END_REF]. Les Nations Unies ont également défini « l'énergie abordable et propre » comme l'un de ses 17 objectifs de développement durable [START_REF]#Envision2030 Goal 7: Affordable and Clean Energy | United Nations Enable[END_REF]. Il vise à accroître l'efficacité énergétique mondiale, la part des énergies renouvelables et à moderniser l'infrastructure énergétique afin que l'accès universel à des services énergétiques abordables, fiables et modernes puisse être assuré d'ici 2030.

Dans le cadre de l'effort mondial et en accord avec les objectifs de l'Accord de Paris [START_REF]Paris Agreement[END_REF], l'Union Européenne (UE) a récemment (2019) adopté un « Paquet Énergie propre pour tous les Européens » [START_REF]Clean energy for all Europeans package[END_REF], initialement proposé en 2016. Avec la vision d'un système énergétique neutre en carbone d'ici 2050, le paquet fixe des objectifs intermédiaires pour l'année 2030 : une réduction de 40% des émissions de gaz à effet de serre (par rapport au niveau de 1990), une augmentation de 32% de la part des sources renouvelables dans le mix énergétique et une efficacité énergétique de 32,5% par rapport à un scénario de statu quo. Le cadre politique ainsi actualisé aide les États membres à renforcer leurs mécanismes de soutien conformément aux nouveaux objectifs. Les divers types de mécanismes de soutien qui ont été pratiqués à travers le continent comprennent les certificats d'électricité verte, les certificats d'économie d'énergie blanche, les prix variables de l'électricité à la fois dans le temps et l'emplacement géographique, les tarifs de rachat, les quotas, les incitations, le net-metering, le bonus-malus, la taxe carbone et ainsi que plusieurs autres leviers [START_REF]IEA -International Energy Agency[END_REF]. La Fig. 1-1 montre les progrès réalisés dans l'intégration des énergies renouvelables dans les secteurs de l'électricité, des transports et du chauffage et du froid dans les pays de l'UE28. Conformément aux tendances depuis 2015, les deux objectifs intermédiaires 2020 fixés pour la part des énergies renouvelables dans le secteur des transports et dans le système énergétique global pourraient ne pas être atteints. On peut voir qu'il y a une intégration relativement plus élevée des énergies renouvelables dans le secteur de l'électricité. L'un des facteurs contributifs est le progrès réalisé dans les technologies telles que l'hydroélectricité, l'éolien et le solaire. De plus, bien que les ressources géothermiques et solaires thermiques soient couramment intégrées dans le secteur du chauffage, elles sont également installées dans le système électrique jusqu'à un certain niveau. En revanche, les technologies basées sur les biocarburants peuvent être intégrées dans les trois secteurs. L'électrification des transports et des secteurs du chauffage et du froid est une alternative pour utiliser tous les potentiels renouvelables afin que la transition énergétique puisse atteindre ses objectifs. Les véhicules électriques et les pompes à chaleur (PAC) sont des exemples de technologies pouvant être utilisées pour électrifier les secteurs du transport et du chauffage, respectivement.

Selon les données de l'Agence internationale de l'énergie, 56,55% de l'électricité et 34,1% de la chaleur consommée dans l'UE28 en 2017 ont été générés à partir de sources non fossiles, dont le nucléaire [START_REF]IEA -International Energy Agency[END_REF]. Cela montre que le secteur du chauffage et du froid est en retard sur le secteur de l'électricité en termes de décarbonation. De plus, le secteur du chauffage et du froid dans l'UE consomme environ 50% de l'énergie finale [11]. Le chauffage des locaux et l'eau chaude sanitaire représentent à eux seuls jusqu'à 79% de la consommation finale d'énergie des ménages de l'UE [START_REF]Heating and cooling (facts and figuers)[END_REF]. Cela implique que des potentiels plus élevés d'intégration et de décarbonation des énergies renouvelables sont disponibles dans le secteur du chauffage et du froid.

La Fig. 1-2 montre la part des énergies renouvelables sur le marché du chauffage et du froid dans les pays membres de l'UE28 en 2018. La Suède est le premier pays générant 65,4% de sa chaleur à partir d'énergies renouvelables, suivie de la Lettonie (55,9%) et de la Finlande (54,6%). Ces pays sont connus pour leur part de marché plus élevée dans le réseau de chauffage urbain (jusqu'à 90% en Suède) avec une part plus élevée de cogénération (jusqu'à 73% en Lettonie) [START_REF]District Energy in Sweden[END_REF][START_REF]District Energy in Latvia | Euroheat & Power[END_REF]. Les Pays-Bas, en revanche, affichent la fraction renouvelable la plus faible (6,1%) du fait que 93% de la demande de chaleur du pays est couverte par l'utilisation de chaudières à gaz [START_REF]District Energy in the Netherlands[END_REF]. La France, avec une part renouvelable de 21,8%, se situe au-dessus de la moyenne de l'UE28 (2018) et même au dessus de l'objectif de l'UE pour 2020 (qui est de 20%).

Le Fig. 1-3 compare la France et la Suède en termes de mix total d'approvisionnement en énergie primaire (TPES) par source pour l'année 2016. Avec un TPES non fossile de 72%, la Suède est l'un des principaux pays en termes de système énergétique durable [START_REF]Energy Policies of IEA Countries -Sweden 2019 Review[END_REF]. Du fait de la part la plus élevée de l'énergie nucléaire, 52% du TPES en France est également non fossile [START_REF]Energy Policies of IEA Countries -France 2016 Review[END_REF]. Cependant, seuls 39% et 10% des TPES en Suède et en France, respectivement, proviennent des énergies renouvelables. Cela montre que la Suède est déjà en avance sur l'objectif d'énergies renouvelables 2030 fixé par l'UE tandis que la France doit faire plus d'efforts.

La France et la Suède figurent parmi les principaux pays à secteur électrique à faible émission de carbone. Comme le montre la Fig. 1-4, 92,3% et 98,6% de l'électricité produite en France et en Suède, respectivement, était produite à partir de technologies à faible émission de carbone en 2016. Le système électrique français est dominé par les centrales nucléaires (72,6% comme le montre la Fig. 1234). La part de l'énergie nucléaire en Suède (40,4% comme le montre la Fig. 1234) est la deuxième plus élevée après la France. La Suède est au troisième rang après le Danemark et la Finlande en termes de carburants à base de biomasse et de valorisation énergétique des déchets [START_REF]Energy Policies of IEA Countries -Sweden 2019 Review[END_REF]. La part de l'éolien, de l'hydroélectricité, des déchets et des biocarburants dans le mix de production d'électricité suédois est supérieure aux valeurs correspondantes en France. En revanche, la part du gaz et du charbon dans le mix de production d'électricité est plus importante en France qu'en Suède. La Fig. 1-5 montre les types de combustibles utilisés pour produire de la chaleur en France et en Suède. La Suède est bien connue pour son secteur du chauffage durable. En France, 54,5% seulement de la chaleur est produite à partir de sources non fossiles, comme le montre la Fig. 1-5(a). Le gaz naturel, couvre 36,9% du secteur français du chauffage. En revanche, 86,5% de la chaleur produite en Suède en 2016 provenait de sources sans fossiles (Fig. 12345). Les biocarburants représentent la part la plus élevée, environ 60,4%, dans le secteur du chauffage suédois.

Énoncé du problème

Idéalement, l'intégration de pompes à chaleur (PAC) à grande échelle et d'un stockage thermique pourrait aboutir à un système énergétique plus durable. Cela signifie que, plus d'énergies renouvelables peuvent être intégrées dans le système électrique sans aucune restriction, les demandes de chaleur de pointe peuvent être efficacement déplacées et, par conséquent, un système énergétique plus efficace, flexible et durable peut être établi. Les résultats basés sur des modèles simplifiés, tels que ceux utilisés dans EnergyPlan [START_REF]EnergyPLAN[END_REF], RETScreen [START_REF]Natural Resources Canada[END_REF] et HOMER [START_REF] Homer Energy | HOMER -Hybrid Renewable and Distributed Generation System Design Software[END_REF], confirment ce potentiel. Dans sa « stratégie de chauffage et de refroidissement », l'UE a également fortement recommandé l'utilisation de chaufferies et de production combinée de chaleur et d'électricité (cogénération) à base de biomasse afin de décarboner et d'augmenter l'efficacité de son secteur du chauffage et du froid [START_REF]Heating and cooling[END_REF][START_REF]EU to fight energy waste with the first Heating and Cooling Strategy[END_REF]. Les pays nordiques sont les premiers à suivre les lignes directrices de l'UE et à intégrer des PAC à grande échelle et des stockages thermiques, comme en témoignent la Suède [36 -38], la Finlande [START_REF] Kontu | Introducing modern heat pumps to existing district heating systems -Global lessons from viable decarbonizing of district heating in Finland[END_REF] et le Danemark [START_REF] Lund | Mapping of potential heat sources for heat pumps for district heating in Denmark[END_REF][START_REF] Lund | Socioeconomic potential for introducing large-scale heat pumps in district heating in Denmark[END_REF].

Cependant, le fonctionnement des chaufferies est uniquement géré par des gestionnaires de réseaux de chaleur sans tenir pleinement compte du réseau de distribution d'électricité. Cela conduit à un fonctionnement sous-optimal du système global. Les gestionnaires de réseaux de distribution d'électricité, par exemple, estiment la demande d'électricité sur la base de ses propres données en considérant les PAC implicitement, généralement comme une demande fixe. Un opérateur de chauffage urbain, en revanche, considère le profil de demande de chaleur et toutes les sources de chaleur, y compris les PAC, et décide de la stratégie opérationnelle la plus économique de ces systèmes. Étant donné que les deux gestionnaires de réseaux de distribution agissent indépendamment, toute différence possible entre la stratégie opérationnelle supposée et réelle des PAC peut compromettre l'efficacité du système global [START_REF] Ma | The optimal structure planning and energy management strategies of smart multi energy systems[END_REF]. De plus, l'insertion accrue des PAC a mis en pause un problème de sécurité car elle surcharge de plus en plus le réseau de distribution d'électricité basse tension [START_REF] Kontu | Introducing modern heat pumps to existing district heating systems -Global lessons from viable decarbonizing of district heating in Finland[END_REF][START_REF] Arnaudo | Techno-economic analysis of demand side flexibility to enable the integration of distributed heat pumps within a Swedish neighborhood[END_REF]. En outre, un placement optimal et une répartition économique des technologies de couplage telles que les cogénérations réalisées uniquement du point de vue du réseau de distribution de chaleur pourraient conduire à des solutions sous-optimales du point de vue du réseau de distribution d'électricité. Cela montre à quel point il est difficile de capturer les détails des paramètres opérationnels du système multi-énergies (MES) à l'aide des outils classiques de modélisation et de simulation de mono vecteur énergétique car ils ne sont pas conçus pour de tels réseaux hybrides [START_REF] Widl | Studying the potential of multi-carrier energy distribution grids: A holistic approach[END_REF].

Des outils de modélisation et de simulation intégrés sont nécessaires pour quantifier cette synergie, évaluer la faisabilité technique et trouver les paramètres opérationnels optimaux d'un système multi-énergies. Les outils de simulation MES simplifiés, comme HOMER et EnergyPlan, n'ont pas cette capacité car ils ne prennent pas en compte les paramètres physiques et de fonctionnement des réseaux de distribution. Une revue de la littérature par Allegrini et al. [START_REF] Allegrini | A review of modelling approaches and tools for the simulation of district-scale energy systems[END_REF] montre qu'il n'existe pas d'outils bien développés pouvant prendre en compte les paramètres physiques des réseaux de distribution MES. Les approches de modélisation sur mesure rapportées dans la littérature (voir l'état de l'art pour plus de détails), en revanche, manquent de généralité et de flexibilité pour répondre aux questions de recherche clés suivantes liées aux réseaux couplés de distribution d'électricité et de chaleur.

2. Comment gérer de manière optimale les technologies de couplage sans compromettre les réseaux deux de distribution ?

3. L'emplacement des technologies de couplage pourrait-il affecter les performances du système énergétique global ? Si oui, comment déterminer leur placement optimal ?

L'état de l'art

Dans la transition vers de futurs programmes d'énergie durable et renouvelable, le MES propose une nouvelle approche pour étudier les solutions possibles exploitant la synergie de différents secteurs de l'énergie, pour aller au-delà des limites de chacun des sous-secteurs individuels [START_REF] Lund | Smart energy and smart energy systems[END_REF][START_REF] Mancarella | MES (multi-energy systems): An overview of concepts and evaluation models[END_REF]. Classiquement, les secteurs de l'énergie sont découplés du point de vue de la planification et de l'exploitation alors qu'en réalité, les deux doivent être abordés conjointement [START_REF] Mancarella | MES (multi-energy systems): An overview of concepts and evaluation models[END_REF][START_REF] Wu | Integrated Energy Systems[END_REF]. Les véhicules électriques, la cogénération et les PAC sont des exemples de technologies énergétiques distribuées qui couplent les secteurs de l'électricité, du gaz, des transports, du chauffage et du froid. Lund et al. [START_REF] Lund | 4th Generation District Heating (4GDH): Integrating smart thermal grids into future sustainable energy systems[END_REF] ont également identifié que la capacité des réseaux thermiques à interagir avec les réseaux électriques intelligents est l'une des caractéristiques des futurs réseaux de chauffage urbain de 4 ème génération. Dans cette optique, l'intégration des centrales de cogénération et des centrales thermiques dans le secteur du chauffage augmente dans les pays nordiques [START_REF] Werner | District heating and cooling in Sweden[END_REF][START_REF] Levihn | CHP and heat pumps to balance renewable power production: Lessons from the district heating network in Stockholm[END_REF][START_REF] Shahrokni | Big meter data analysis of the energy efficiency potential in Stockholm's building stock[END_REF][START_REF] Kontu | Introducing modern heat pumps to existing district heating systems -Global lessons from viable decarbonizing of district heating in Finland[END_REF][START_REF] Lund | Mapping of potential heat sources for heat pumps for district heating in Denmark[END_REF][START_REF] Lund | Socioeconomic potential for introducing large-scale heat pumps in district heating in Denmark[END_REF], qui sont des modèles pour des systèmes énergétiques de quartier durables. Une PAC utilise l'électricité pour élever une température basse à une température plus élevée, tandis qu'un cogénération produit à la fois de la chaleur et de l'électricité à partir de sources de carburant, telles que les déchets, la biomasse et le gaz naturel. Intégrer davantage ces technologies vise à renforcer le couplage et la synergie entre les réseaux.

Cependant, afin de quantifier cette synergie et de gérer efficacement les systèmes énergétiques de couplage, des outils de modélisation et de simulation appropriés sont nécessaires. Le fait d'avoir un certain nombre de dispositifs de conversion et de stockage d'énergie avec des caractéristiques différentes qui impliquent plusieurs vecteurs énergétiques et de nombreux paramètres rend le processus de modélisation et de simulation d'un MES très complexe et exigeant [START_REF] Ma | The optimal structure planning and energy management strategies of smart multi energy systems[END_REF]. La prise en compte des paramètres de fonctionnement importants d'un MES à l'aide des outils conventionnels est lourde car ces outils ne sont pas conçus pour des systèmes de distribution d'énergie multi-vecteurs énergétiques [START_REF] Widl | Studying the potential of multi-carrier energy distribution grids: A holistic approach[END_REF]. Bien qu'il n'existe pas encore d'outil logiciel disponible sur le marché capable de capturer l'intégralité des paramètres physiques et de fonctionnement des réseaux de distribution multi-opérateurs [START_REF] Allegrini | A review of modelling approaches and tools for the simulation of district-scale energy systems[END_REF][START_REF] Van Beuzekom | A review of multi-energy system planning and optimisation tools for sustainable urban development[END_REF], des progrès ont été réalisés récemment au niveau de la recherche.

L'ensemble des travaux présentés dans la littérature peut être largement regroupé en deux grands groupes selon que les détails des paramètres du réseau de distribution sont pris en compte ou non (voir Fig.

Méthodologie

La structure de l'ensemble du modèle développé dans cette thèse est illustrée à la Fig.

Résultats et discussions sur l'intégration des études de flux de charge

Les études de cas de démonstration sont soigneusement conçues pour illustrer la capacité de l'outil développé à combler les lacunes scientifiques discutées précédemment. Tout d'abord, l'outil est utilisé pour mener une étude statique d'une heure sur le flux de charge sur les réseaux d'électricité et de chauffage couplés à plusieurs centrales de production combinée de chaleur et d'électricité (CHP). Quatre études de cas avec des topologies à la fois radiales et maillées sont considérées comme le montre la Fig. 5-2.

Un des cas (Case 5-A sur la Fig. 5-2a) est considéré pour comparaison avec les travaux antérieurs rapportés dans la littérature [START_REF] Shabanpour-Haghighi | An Integrated Steady-State Operation Assessment of Electrical, Natural Gas, and District Heating Networks[END_REF]. Les résultats de l'étude de cas comparative sont présentés dans le Les profils de température et de pression hydraulique après la solution de flux de charge sont représentés sur la Fig. 5-4, tandis que les profils de tension sont illustrés sur la Fig. 5-5.

Le potentiel de l'outil développé pour effectuer une simulation pseudo-dynamique est démontrée en considérant une étude de cas à six hubs (comme le montre la Fig. 56). L'année est représentée par huit jours, 6 février, 23 mars, 6 mai, 23 juin, 7 août, 23 septembre, 7 novembre et 23 décembre, donnant un total de 192 heures de simulation. Ces profils sont pris uniquement à titre d'illustration pour démontrer la capacité de l'outil proposé à gérer la simulation pseudo-dynamique d'un MES fortement couplé. Le Table 5-10 résume le profil de demande annuel pour chaque vecteur pour l'ensemble du réseau. L'énergie solaire photovoltaïque, l'incinération des déchets par cogénération, le gaz pour la chaudière et les PACs sont considérés pour satisfaire la demande de chaleur et d'électricité. La variation annuelle des paramètres opérationnels, notamment l'amplitude de la tension, l'angle de tension, le débit massique, les températures d'alimentation et de retour des hubs, les températures d'entrée et de sortie de l'eau dans les conduites aller et retour, les pertes de puissance sur les réseaux d'électricité et de chauffage, la consommation de combustible, la quantité d'énergie importée et exportée sont simulées en fonction du temps.

Les résultats de l'ensemble, une heure et plusieurs heures (pseudodynamique) simulation, ont démontré la flexibilité du modèle proposé. La capacité du modèle à capter tous les paramètres de fonctionnement importants des réseaux de distribution de chaleur et d'électricité est démontrée. Les inversions de débit sont automatiquement traitées par l'algorithme. Le modèle proposé considère également la variation du coefficient de transfert de chaleur global en conduite en fonction du débit massique et de la température du sol.

Son évolutivité par sa capacité à intégrer différentes technologies de générations et de couplage distribuées via les matrices de couplage est également illustrée.

Les simulations ont également été exécutées avec des conditions initiales différentes et ont convergé vers la même solution, mais avec un nombre différent d'itérations. Il est également démontré que le modèle proposé est plus précis et réaliste par rapport à ce qui a été précédemment rapporté dans la littérature.

Une comparaison entre les réseaux radiaux et maillés (pour des simulations horaires) a révélé la possibilité de pertes de distribution plus élevées dans le cas maillé. D'autre part, l'étude de cas de simulation pseudo-dynamique montre que la façon dont la cogénération et les PACs fonctionnent affecte aussi bien les réseaux d'électricité et de distribution de chaleur. Le dimensionnement des capacités de production des technologies de couplage en fonction d'un seul des deux types de profil de demande pourrait conduire à une capacité de production sous-dimensionnée ou surdimensionnée pour l'autre réseau. Ce travail présente également des pertes de distribution et le coût total d'exploitation qui varient en fonction de l'emplacement des pompes à chaleur. Cela montre que la prise en compte simultanée des deux réseaux est cruciale pour optimiser les technologies de couplage telles que les HP.

Résultats et discussions sur le flux de puissance optimal intégré

Une technique d'optimisation des essaims de particules (particle swarm optimisation -PSO) est intégrée au problème de flux de charge développé, et est appliquée pour trouver le placement optimal des technologies de couplage, telles que les HP et les chaudières. La fonction objective correspond à la minimisation du coût d'exploitation. Un algorithme PSO imbriqué (illustré à la Fig. 6-1) est ensuite appliqué à deux cas d'étude correspondants à deux tailles de structures de réseaux (Fig. 6-2 et Fig. 6-6, respectivement).

L'emplacement optimal et le dimensionnement des PACs et des chaudières pour les petites et les grandes études de cas sont résumés dans le Table 6-2 et le Table 6-4, respectivement. Le débit massique injecté dans le réseau à partir de chaque hub et les températures aller et retour correspondantes pour le cas de la petite structure de réseaux sont illustrés à la Fig. 6-3. Des résultats similaires pour l'étude de la structure de réseaux plus grande sont présentés sur la Fig. 67. Le débit massique dans chaque conduite et les tensions à chaque hub sont présentés sur les Fig. 6 L'étude d'un cas supplémentaire (indiqué sur la Fig. 6-12) est considéré pour en vue d'une optimisation de la structure de distribution par identification des branches des réseaux les moins performantes. Une approche étape par étape est présentée dans ce sens Fig. 6-11 et la comparaison des pertes de distribution à différentes étapes est présentée sur la Fig. 67891011121314.

Les résultats montrent que l'optimisation de l'emplacement et de la taille des PAC uniquement du point de vue du réseau de chauffage pourrait réduire les pertes de chaleur, mais parfois elle pourrait surcharger et augmenter les pertes dans le réseau de distribution d'électricité, en particulier lorsque les HP sont installés sur des réseaux de distribution basse tension. Cela pourrait également conduire à des profils de tension inacceptables. D'autre part, une optimisation intégrée qui prend en compte les contraintes et les pertes de distribution dans les deux réseaux se traduit par un coût d'exploitation inférieur pour l'ensemble du système avec la possibilité d'un meilleur profil de tension dans le réseau de distribution d'électricité.

On constate que les pertes de charge par frottement sont tout aussi importantes que les pertes de charge aux postes consommateurs pour les réseaux plus importants. L'optimisation thermo-économique a en outre révélé que l'abaissement de la température d'alimentation n'entraîne pas nécessairement une baisse des coûts d'exploitation. Le coût de pompage peut parfois contrebalancer les avantages de températures d'alimentation plus basses. Il convient donc de prendre en compte le coût total du pompage dans les optimisations thermoéconomiques des grands réseaux. En revanche, l'abaissement de la température de retour des consommateurs (côté primaire) s'avère toujours économique.

Résultats et discussions sur l'intégration optimale du stockage thermique

Un modèle de réservoir d'eau de stockage thermique sensible stratifié est intégré au modèle d'optimisation. Un nouvel algorithme de contrôle prédictif basé sur un modèle PSO (MPCillustré sur la Fig. 7-3a) est mis en oeuvre pour analyser la synergie entre les réseaux d'électricité et de chauffage en présence d'énergies renouvelables intermittentes et d'un signal de prix de l'électricité variable. Une étude de cas illustrée à la Fig. 7-4 est utilisée avec les données d'entrée résumées aux Fig. 7-5 et Fig. 7-6.

La Fig. 7-7 montre le bilan énergétique de la chaleur, tandis que la Fig. 7-8 illustre l'état de charge et les paramètres de fonctionnement du stockage thermique. Les températures aller et de des hubs sélectionnés sont indiquées sur la Fig. 789. Tandis que les débits des conduites sont illustrés à la Fig. 7-10. De même, le bilan électrique actif, le bilan électrique réactif, les profils de tension et les courants de ligne sont représentés sur les Fig. 7-11, Fig. 7-12, Fig. 7-13 et Fig. 7-14, respectivement. Le Table 7-2 compare le système avec et sans le stockage thermique en utilisant les paramètres clés du système énergétique global. Les deux systèmes ont utilisé une quantité importante de l'électricité excédentaire générée par la centrale éolienne.

En accord avec [START_REF] Vandermeulen | Controlling district heating and cooling networks to unlock flexibility: A review[END_REF], la perte de distribution dans le réseau de chauffage avec stockage thermique est supérieure à celle du système sans stockage (voir Table 7-2). La raison en est l'exigence d'une température plus élevée pour charger le stockage. La consommation électrique plus élevée des pompes de circulation dans le système sans stockage indique également la présence d'un débit plus élevé dans les conduites en raison de températures d'alimentation plus basses. Cependant, il n'y a pas de différence significative sur les pertes de distribution d'électricité.

Le système avec le stockage thermique évitait totalement l'utilisation de la chaudière à gaz en déplaçant les demandes de chaleur de pointe dans les heures bon marché du signal de prix de l'électricité. Cela a réduit la consommation de gaz dans le système et le coût fonctionnel global.

Conclusions

Cette thèse présente une nouvelle approche en modélisation, simulation et optimisation de systèmes multi-énergies avec une prise en compte complète des paramètres de fonctionnement des réseaux de distribution multi-vecteurs. Un cadre général de modélisation du problème de flux de charge en régime permanent des systèmes énergétiques multi-vecteurs est formulé sur la base du concept de hub énergétique étendu. L'approche de modélisation est ensuite illustrée en développant un outil qui peut être utilisé pour mener des études intégrées de flux de charge sur des réseaux d'électricité et de chauffage fortement couplés.

L'algorithme, développé sur MATLAB, se compose d'un modèle hydraulique plus flexible et d'un modèle thermique plus réaliste qui le rend suffisamment général pour que chaque fois qu'il est nécessaire de simuler différentes configurations de réseau, l'utilisateur n'a qu'à alimenter les données du réseau. Le système d'équations et la matrice Jacobienne nécessaire pour résoudre le problème de flux de charge intégré sont automatiquement configurés dans cet outil qui permet de surmonter les limites dans les modèles présentés dans la littérature. Une étude comparative a également montré que le modèle proposé dans cette thèse a une meilleure précision.

Tout changement possible de la direction du débit massique est automatiquement traité dans l'algorithme, comme le confirment les différentes études de cas. L'outil considère également la variation du coefficient de transfert de chaleur global d'un tuyau en fonction du débit massique, ce qui le rend plus réaliste. L'outil développé est également très flexible de telle sorte que différentes générations réparties peuvent être facilement intégrées en tant que générations locales.

La capacité de l'outil est également démontrée dans la gestion des flux d'énergie bidirectionnels à travers les branches et entre les concentrateurs et le réseau. Un hub énergétique donné peut être un « prosommateur » (producteur à certains instants et consommateur à d'autres) et non plus seulement consommateur. Différentes technologies de couplage peuvent également être intégrées dans l'outil via les matrices de couplage. Les cogénérations, les HPS, les chaudières, les centrales éoliennes, les stockages, les panneaux solaires photovoltaïques et thermiques sont pris en compte dans différentes études de cas en tenant pleinement compte des paramètres de fonctionnement des réseaux de chauffage et de distribution d'électricité.

L'outil a été testé pour différentes topologies de réseaux constituées de réseaux d'électricité et de chauffage fortement couplés. Les résultats trouvés sont cohérents avec la théorie. Les simulations ont également été effectuées dans des conditions initiales différentes et ont convergé vers la même solution, mais avec un nombre d'itérations différent. Les résultats de l'étude de flux de charge intégrée donnent les détails des paramètres du réseau tels que la tension, le débit massique, la température, la pression, les pertes thermiques, les pertes électriques, les générations aux noeuds d'équilibrage et la quantité de flux d'énergie à travers chaque conduite et ligne de transmission. Les études de cas présentées montrent que les pertes et les paramètres du réseau sont fortement affectés par la topologie des réseaux.

Des études de flux de charge sur plusieurs heures ont également été menées sous la forme d'une simulation pseudo-dynamique sur des réseaux électriques et de chauffage couplés comprenant une cogénération, un solaire PV, une chaudière à gaz et des PACs. Il a été observé que le fonctionnement de la cogénération et des pompes à chaleur affecte les deux réseaux. En outre, le dimensionnement des technologies de production et de couplage sur la base d'un seul type de profil de demande pourrait conduire à une capacité de production sousdimensionnée ou surdimensionnée pour l'autre réseau. Des études comparatives ont également montré que l'emplacement des HP affecte les pertes de distribution et les paramètres de fonctionnement dans les deux réseaux de distribution.

Pour une répartition économique des unités de production et des technologies de couplage, en tenant pleinement compte des paramètres des réseaux de distribution, un algorithme d'optimisation méta-heuristique a été développé sur la base de la technique PSO. Il est appliqué pour trouver l'emplacement et les tailles optimales des technologies de génération et de couplage distribuées. Une étude comparative est faite entre la gestion uniquement du point de vue du réseau de chauffage et la gestion en considérant les deux réseaux. Les résultats montrent que l'optimisation de l'emplacement et de la taille des HP uniquement du point de vue du réseau de chauffage pourrait réduire les pertes de chaleur. Mais, il est démontré qu'une telle stratégie de gestion peut entraîner une surcharge et davantage de pertes dans le réseau de distribution d'électricité. Cela pourrait également conduire à des profils de tension inacceptables. Le problème s'aggrave, en particulier lorsque de grandes PACs sont installées dans un réseau de distribution basse tension. D'autre part, une optimisation intégrée qui prend en compte les contraintes et les pertes de distribution dans les deux réseaux a entraîné une baisse des coûts d'exploitation pour l'ensemble du système et un meilleur profil de tension dans le réseau électrique.

Les pompes de circulation sont utilisées pour surmonter deux types de pertes de charge. Le premier correspond aux pertes par frottement dans les canalisations tandis que le second est lié à la différence de pression minimale requise aux postes (côté primaire). Leur ampleur relative dépend de la taille et de la stratégie d'exploitation du réseau de chauffage. On constate que la consommation électrique des pompes de circulation n'est généralement pas significative par rapport à la demande totale d'électricité, pour les réseaux étudiés. Cependant, comme le coût total du pompage de l'énergie peut concurrencer le coût de la perte de distribution, surtout s'il existe des sources de chaleur très bon marché, le coût du pompage doit être pris en compte dans les optimisations thermo-économiques.

Une approche étape par étape pour identifier et isoler les branches avec perte d'un réseau de chaleur maillé a été réalisée à l'aide d'un flux de puissance optimal intégré et de procédures d'optimisation thermo-économique. Une telle approche permet de faire fonctionner le réseau de chaleur efficacement en dehors des heures de pointe. Les résultats de l'optimisation thermo-économique ont en outre révélé que l'abaissement de la température d'alimentation n'entraîne pas nécessairement un coût d'exploitation inférieur. Le coût de pompage peut parfois contrebalancer les avantages de températures d'alimentation plus basses. Cependant, l'abaissement de la température de retour des consommateurs (côté primaire) s'avère toujours économique.

La synergie des réseaux de chauffage et d'électricité couplés peut être davantage exploitée à l'aide de réservoirs de stockage thermique. Les stockages thermiques peuvent augmenter la flexibilité du réseau de chaleur en réduisant les pics de demande. Ils augmentent également la pénétration des énergies renouvelables dans les réseaux de distribution d'électricité et de chaleur. Un système énergétique comprenant une cogénération, une PAC, une chaudière à gaz, une éolienne et un stockage thermique a été étudié. Un signal de prix variable pour l'électricité du réseau est également pris en compte. Un MPC basé sur PSO est mis en oeuvre pour optimiser le système sur un horizon de 24 heures. Il est démontré qu'une part importante de la production éolienne peut être utilisée dans le système en couplant les réseaux d'électricité et de chauffage. Le système avec le stockage thermique évitait totalement l'utilisation de la chaudière à gaz en déplaçant les demandes de chaleur de pointe dans les heures bon marché du signal de prix de l'électricité. Cela a réduit la consommation de gaz dans le système et le coût fonctionnel global.

Il est démontré (à travers des études de cas) que les méthodologies de simulation et d'optimisation basées sur la notion de hubs énergétiques étendus, proposées dans cette thèse sont évolutives, flexibles et très efficaces pour considérer les paramètres de fonctionnement importants et quantifier la synergie qui peut être obtenue à partir des réseaux intégrés d'électricité et de chauffage. La flexibilité supplémentaire qui peut être obtenue en intégrant le stockage thermique est également explorée en tenant compte de toutes les contraintes techniques.

Mots clés :

Réseaux de chauffage urbain; Réseau de distribution d'électricité; Hub énergétique étendu; Étude intégrée du flux de charge; Flux de puissances intégrées optimaux; Réseaux énergétiques multi-vecteurs; Système multi-énergies; Placement optimal; Stockage thermique.

xx | Doctoral thesis by Getnet Tadesse Ayele

Sammanfattning

Ny litteratur visar att det finns en betydande potential för avkolning och effektivitetsförbättring som kan uppnås genom synergi från multi-energisystem (MES). Kopplingsteknik, såsom kraftvärmeanläggningar, värmepumpar och värmelager rekommenderas allmänt som medel för att låsa upp ytterligare flexibilitet och öka penetrationen av förnybara energikällor inom värme-och elsektorn. Mot bakgrund av detta ökar storleken och antalet kopplingstekniker, såsom kraftvärmeverk och värmepumpar (HP), som installeras i värmefördelningsnätet. Eftersom dessa tekniker uteslutande hanteras av fjärrvärmenätoperatörerna blir deras drift ibland suboptimal ur elnätets synvinkel, och de (särskilt stora HP) kan orsaka överbelastning av distributionsnätet för lågspänning.

Integrerade simulerings-och optimeringsmodeller krävs för att effektivt utnyttja synergierna utan att kompromissa med MES: s ingående distributionsnät. Sådana modeller är ännu inte väl utvecklade. De konventionella simuleringsverktygen för enenergibärare kan inte heller fånga viktiga driftsparametrar i distributionsnätet med flera bärare.

En ny metod för simulering och optimering av MES har utvecklats i denna avhandling baserad på ett EEH-synsätt (Extended Energy Hub). Det allmänna ramverket utvecklas först i modulär form så att det enkelt kan anpassas för alla typer av multibärare energinät. Ramverket används sedan för att utveckla detaljerna i en integrerad belastningsflödesmodell som styr kopplade värme-och eldistributionsnät. Olika fallflödesstudier med radiella och meshade topologier övervägs för demonstration och numerisk validering av den föreslagna modellen.

Lastflödesmodellen kombineras vidare med en algoritm för optimering av partikelvärm för att genomföra integrerade optimala effektflödesstudier. Dess bidrag till teknikens ståndpunkt demonstreras genom att studera den optimala placeringen av kopplingstekniker, såsom HP och pannor i kopplade värme-och eldistributionsnät. Modellens kapacitet illustreras ytterligare genom att utnyttja synergierna med hjälp av HP tillsammans med termisk lagring i närvaro av intermittenta förnybara energikällor och variabel elprissignal.

Det visas att de EEH-baserade simulerings-och optimeringsmetoderna som föreslås i denna avhandling är mycket effektiva, flexibla och lätt skalbara för att fånga de viktigaste driftsparametrarna för integrerade el-och fjärrvärmenät. Modellerna kan användas som en plattform för vidare studier om integration av smarta nät och smarta termiska nätverk. The research focuses on optimal integration of electricity and heating networks. It was conducted in two (home and host) institutes. The Department of Energy Systems and Environment (DSEE) at IMT Atlantique, Nantes, France (which is part of the UMR CNRS 6144 GEPEA) served as a home institute, and it is where slightly more than half of the research period was spent. The remaining time was used at the host institute in the Division of Heat and Power (EKV), Department of Energy Technology (EGI), School of Industrial Engineering and Management (ITM), KTH Royal Institute of Technology, Stockholm, Sweden.

Nyckelord

As a wrap-up to the research conducted, this thesis presents new modelling, simulation and optimisation approaches that can be used to exploit the synergies from coupled heating and electricity networks without compromising the respective distribution networks.

Nantes, France (Dec. 2020)

Getnet Tadesse Ayele Thermal conductivity of the carrier pipe (W.m -1 .K -1 ) 𝑘 3

I
Thermal conductivity of the insulating material of the pipe (W.m -1 .K -1 ) 𝑘 [START_REF]IEA -International Energy Agency[END_REF] Thermal conductivity of the outer jacket of the pipe (W.m -1 .K -1 ) 𝑘 𝑠 Thermal conductivity of the soil (W.m -1 .K -1 ) 𝑘 𝑤 Thermal conductivity of the water (W.m - approximations are acceptable only if the distribution network is small or if the reactive power flow in the system is not significant.

List of Tables

List of Figures

Li et al. [START_REF] Li | Optimal operation of the integrated electrical and heating systems to accommodate the intermittent renewable sources[END_REF] studied the electricity and heating networks that are coupled with a HP. AC power flow model is applied for electricity network while mass balance and temperature drops are used to model the heating network. The possibility of meshed heating network topologies is not considered. Also, a single coupling device is assumed. Awad et al. [START_REF] Awad | Integrated optimal power flow for electric power and heat in a MicroGrid[END_REF] and Liu et al. [START_REF] Liu | Combined analysis of electricity and heat networks[END_REF] followed a similar approach, but node-loop equations with pseudo-dynamic temperature drop equations are used to model the heating network considering meshed loops. However, the overall heat transfer coefficient is assumed to be constant for all mass flows, which is not true in reality. Furthermore, the node-loop equations used in the hydraulic models need pseudo-loop paths (the number of which depends on the network topology) in addition to the physical loops in order to have equal number of equations as the number of unknowns. As identification of pseudo-loops is not straight forward, it is difficult to develop a general algorithm for the models that are based on node-loop equations [START_REF] Boulos | Comprehensive Water Distribution Systems Analysis Handbook for Engineers and Planners[END_REF]. In their case study, Liu et al. [START_REF] Liu | Combined analysis of electricity and heat networks[END_REF] considered heating and electricity networks that are coupled through three CHP plants. All electrical demands and all CHP plants are assumed to be at unity power factor which is not always the case. There is also a requirement of reactive power to keep the voltage magnitude at the slack and PV hubs to a specified magnitude. However, no reactive power generation is considered in their case study [START_REF] Liu | Combined analysis of electricity and heat networks[END_REF].

Liu and Mancarella [START_REF] Liu | Modelling, assessment and Sankey diagrams of integrated electricity-heat-gas networks in multi-vector district energy systems[END_REF] and Shabanpour-Haghighi and Seifi [START_REF] Shabanpour-Haghighi | An Integrated Steady-State Operation Assessment of Electrical, Natural Gas, and District Heating Networks[END_REF] extended the load flow model presented in [START_REF] Liu | Combined analysis of electricity and heat networks[END_REF] into an energy system consisting of gas, electricity and DHN. The application of such models to study network congestion and uncertainties are presented in [START_REF] Martinez Cesena | Energy Systems Integration in Smart Districts: Robust optimisation of Multi-Energy Flows in Integrated Electricity, Heat and Gas Networks[END_REF] and [START_REF] Khorsand | Probabilistic energy flow for multi-carrier energy systems[END_REF], respectively. Both active and reactive power flows are considered in their case studies. However, lack of generality on the hydraulic models remained unsolved due to the use of the node-loop equations. They also assumed constant value of the overall heat transfer coefficients in their thermal models.

In a load flow study of a MES, there are a number of equations for each network to be solved iteratively. Due to the coupling devices, some of the equations in one type of energy network are coupled with some of the equations in the other type of energy network. To handle this coupling issue, two types of iterative strategies are reported in literature: a decomposed approach and an integrated approach. In the first method, the overall system is decomposed into independent network models and one network model is solved at a time based on the previous values of the parameters in the other networks [START_REF] Rees | Steady state flow analysis for integrated urban heat and power distribution networks[END_REF][START_REF] Moeini-Aghtaie | A Decomposed Solution to Multiple-Energy Carriers Optimal Power Flow[END_REF]. The integrated approach, on the other hand, considers the coupling between different networks explicitly and solves the overall system of equations as a single load flow problem as described in [START_REF] Shabanpour-Haghighi | An Integrated Steady-State Operation Assessment of Electrical, Natural Gas, and District Heating Networks[END_REF][START_REF] Batić | Combined energy hub optimisation and demand side management for buildings[END_REF]. A comparative study of the two approaches presented in reference [START_REF] Liu | Combined analysis of electricity and heat networks[END_REF] shows that the integrated load flow approach is faster than the decomposed approach and the difference in computational time becomes significant as the network size gets larger.

There are well-developed simulation and optimisation tools, on the other hand, capable of handling only one energy carrier. Quite few researchers combined two and more of them in a co-simulation environment so that distribution networks' operating parameters can be investigated. Their application for integrated optimisation is yet to be discovered.

Various tailor-made MES models are reported in literature with and without the consideration of distribution networks. Most of them are case-specific while quite a few followed a modular approach. The modular tailor-made models applied the energy hub concept together with network power flows. However, the unidirectional power flow assumption at the energy hubs and isothermal hydraulic models limits their application scope to a MES consisting of only electricity and gas networks. The models cannot handle prosumers either.

The energy hub concept shall, therefore, be extended (redefined) to add a bidirectional power flow flexibility and increase its capacity to handle heat distribution networks. Also, the thermo-hydraulic models discussed in the tailor-made approaches need to be general enough to include meshed topologies. In addition, the models are preferred to be more flexible so that any change in the topology or configuration of the heating network can be automatically detected in the computer-based algorithm. Furthermore, the modelling approach should treat the MES as a single problem so that multi-carrier power flows can be simulated simultaneously. This, in turn, makes it easier to apply an optimisation technique for the whole system. The scalability and modularity of the methodology to encompass any additional energy carrier or coupling technology is also crucial. The following chapter introduces a novel modelling framework that can address all these issues. matrix representation of the relationship between the demands (𝑚 types), the net power inputs into the coupling system (𝑛 types) and the rate of change of energy in the storage units (𝑛 types) [START_REF] Geidl | Integrated Modeling and Optimisation of Multi-Carrier Energy Systems[END_REF].

[

𝐿 𝑎 ⋮ 𝐿 𝛿 ⋮ 𝐿 𝑚 ] = [ 𝐶 𝑎𝑎 ⋯ 𝐶 𝑎𝛾 ⋯ 𝐶 𝑎𝑛 𝑍 𝑎𝑎 ⋯ 𝑍 𝑎𝛾 ⋯ 𝑍 𝑎𝑛 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 𝐶 𝛿𝑎 ⋯ 𝐶 𝛿𝛾 ⋯ 𝐶 𝛿𝑚 𝑍 𝛿𝑎 ⋯ 𝑍 𝛿𝛾 ⋯ 𝑍 𝛿𝑛 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 𝐶 𝑚𝑎 ⋯ 𝐶 𝑚𝛾 ⋯ 𝐶 𝑚𝑚 𝑍 𝑚𝑎 ⋯ 𝑍 𝑚𝛾 ⋯ 𝑍 𝑚𝑛] [ 𝑃 𝑎_𝑖𝑛 ⋮ 𝑃 𝛾_𝑖𝑛 ⋮ 𝑃 𝑛_𝑖𝑛 𝐸 ̇𝑎 ⋮ 𝐸 ̇𝛾 ⋮ 𝐸 ̇𝑛 ] (3-1)
where 𝐿 𝛿 is the local demand of type 𝛿; 𝑃 𝛾_𝑖𝑛 is the net power input into the coupling system of type 𝛾 energy carrier; 𝐸 ̇𝛾 is the time rate of change of stored energy of type 𝛾; 𝐶 𝛿𝛾 is an element in the coupling matrix showing how much of the input power to the coupling system of type 𝛾 is converted into the local demand type of 𝛿; 𝑍 𝛿𝛾 is an element of the coupling matrix indicating how much of the power from a storage type 𝛾 is contributing to meet the local demand of type 𝛿; 𝛿 ranges from a to m while 𝛾 ranges from a to n.

By referring Fig. 3-2, it can be seen that each of the net input powers into the coupling system, (𝑃 𝛾_𝑖𝑛 ) terms in equation (3-1), is equal to the corresponding local generation (𝑃 𝛾𝑔 ) minus the injection into the network (𝑃 𝛾 ) of the corresponding energy carrier (see equation (3-2)). By substituting equation (3-2) into equation (3-1), the extended version of the coupling equation for a general network of MES can be formulated as shown in equation . As the power injections into the network (which are now part of the coupling equation) are functions of the network parameters at all hubs, they act as a link between the demands and generations at different hubs.

𝑃 𝛾_𝑖𝑛 = 𝑃 𝛾𝑔 -𝑃 𝛾 (3-2) [ 𝐿 𝑎 ⋮ 𝐿 𝛿 ⋮ 𝐿 𝑚 ] = [ 𝐶 𝑎𝑎 ⋯ 𝐶 𝑎𝛾 ⋯ 𝐶 𝑎𝑛 𝑍 𝑎𝑎 ⋯ 𝑍 𝑎𝛾 ⋯ 𝑍 𝑎𝑛 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 𝐶 𝛿𝑎 ⋯ 𝐶 𝛿𝛾 ⋯ 𝐶 𝛿𝑚 𝑍 𝛿𝑎 ⋯ 𝑍 𝛿𝛾 ⋯ 𝑍 𝛿𝑛 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 𝐶 𝑚𝑎 ⋯ 𝐶 𝑚𝛾 ⋯ 𝐶 𝑚𝑚 𝑍 𝑚𝑎 ⋯ 𝑍 𝑚𝛾 ⋯ 𝑍 𝑚𝑛] [ 𝑃 𝑎𝑔 -𝑃 𝑎 ⋮ 𝑃 𝛾𝑔 -𝑃 𝛾 ⋮ 𝑃 𝑛𝑔 -𝑃 𝑛 𝐸 ̇𝑎 ⋮ 𝐸 ̇𝛾 ⋮ 𝐸 ̇𝑛 ] (3-3)
The different numerical methods that can be used to solve hydraulic head and mass flow rates in water supply networks are presented in [START_REF] Boulos | Comprehensive Water Distribution Systems Analysis Handbook for Engineers and Planners[END_REF][START_REF] Larock | Hydraulics of Pipeline Systems[END_REF], while their applications in DHNs are discussed in [START_REF] Pirouti | Modelling and analysis of a district heating network[END_REF]. The methods range from a classical Hardy-Cross loop-correction method to a computer based Newton-Raphson method. In any of the methods, continuity of flow (conservation of mass at each node) and conservation of energy (pressure drop equations) are the principles that should be guaranteed [START_REF] Nielsen | Methods for analyzing pipe networks[END_REF]. Boulos et al. [START_REF] Boulos | Comprehensive Water Distribution Systems Analysis Handbook for Engineers and Planners[END_REF] discussed four types of methods that can be used to combine and solve the two principles. These are loop equations, node-loop equations, node equations and pipe equations.

The loop equations are formulated using a non-linear relationship between pressure drop and mass flow rate across a closed loop. Although this method could give the least possible number of equations, it requires initial values that satisfy the conservation of mass throughout the network. Moreover, additional pseudo-loops (the number of which depends on the topology) are usually needed to balance the number of equations and the number of unknowns, which makes it difficult to use them on computer-based algorithms.

The node-loop equations method, on the other hand, uses mass conservation equations at each node and loop equations for closed loops. Additional pseudo-loops are defined, depending on the topology. This method doesn't need precise initial guess of mass flow rates, but the necessity of pseudo-loops is still a drawback. Due to this inflexibility, the node-loop equations are used to model only specific case studies of DHN as presented in [START_REF] Awad | Integrated optimal power flow for electric power and heat in a MicroGrid[END_REF][START_REF] Liu | Modelling, assessment and Sankey diagrams of integrated electricity-heat-gas networks in multi-vector district energy systems[END_REF][START_REF] Liu | Combined analysis of electricity and heat networks[END_REF][START_REF] Shabanpour-Haghighi | An Integrated Steady-State Operation Assessment of Electrical, Natural Gas, and District Heating Networks[END_REF][START_REF] Rees | Steady state flow analysis for integrated urban heat and power distribution networks[END_REF][START_REF] Pirouti | Modelling and analysis of a district heating network[END_REF][START_REF] Fedorov | Parallel implementation of a steady state thermal and hydraulic analysis of heat exchanger[END_REF].

In contrast to the previous two methods, in which every equation is expressed as a function of pipe mass flow rate, the node equations method defines the conservation of mass at each node as a function of nodal pressures. Hence, it requires only a number of equations equal to the number of nodes. This is analogous to the electricity network equations which are based on nodal voltage profiles. However, since the mass flow rates are the coupling parameters between the hydraulic and thermal models, expressing them as a highly non-linear function of pressure head makes the overall system highly complex and difficult to handle.

Unlike the other three methods, the pipe equation uses conservation of mass at each node and pressure drop equations for each pipe. Hence, both parameters (pressure and mass flow rate) are involved in this method with the cost of additional equations. However, the equations involved are relatively easier to compute. As a result, hydraulic models based on pipe equations result in a simple to handle model with equal (even better) convergence rate when compared to the other methods [START_REF] Boulos | Comprehensive Water Distribution Systems Analysis Handbook for Engineers and Planners[END_REF]. Hence, this method is used to represent the hydraulic part of the heating network.

Assuming incompressible and constant water properties throughout the network, the conservation of mass and pressure drop equations can be written as shown in equations (4-8) and (4-9) respectively. These equations are valid for both supply and return pipe networks of the heating network. It is assumed that each of the pipes in the network are circular and straight with uniform cross-section throughout their length. It is reported that the Darcy-Weisbach equation for the pressure resistance coefficient, 𝑘 𝑖𝑗 , is applicable for a wide range of mass flow rate [START_REF] Larock | Hydraulics of Pipeline Systems[END_REF]. Accordingly, 𝛽 (the exponential constant) is taken to be 2 and 𝑘 𝑖𝑗 is given by equation (4-10). where 𝐻 𝑗 and 𝐻 𝑖 are the hydraulic heads (due to pressure and elevation) at nodes 𝑗 and 𝑖 respectively while 𝑚̇𝑗 𝑖 is the mass flow rate flowing from node j to node 𝑖; 𝐿 is the length of the pipe with internal diameter 𝐷 while 𝜌 is the density of water and f is the friction factor. Equation (4-11) can be used to find the value of the friction factor (𝑓) for laminar flow as a function of Reynolds number (𝑅 𝑒 ), absolute surface roughness (𝜀) and pipe diameter (𝐷). For turbulent flow (𝑅 𝑒 > 4000), the exact value of the friction factor is determined using the Colebrook equation as indicated in equation ( 4 conductive shape factor (𝑆). For very long (relative to the diameter) underground pipes, the conductive shape factor per unit length is given by equation (4-17) [START_REF] Bergman | Fundamentals of heat and mass transfer[END_REF]. If the soil temperature around the pipe is known, then equation (4-18) can be used to calculate the thermal resistance between the water inside the pipe and the outer surface of the casing jacket.

As both of them are approximations of the heat transfer between a pipe and a semi-infinite sphere (earth), they have their own strength and limitation. Equation (4-16) is more accurate with the assumption of constant surface temperature, but does not take the earth's internal heat source into account. Equation [START_REF]IEA -International Energy Agency[END_REF][START_REF] Mohammadi | Energy hub: From a model to a concept -A review[END_REF][START_REF]Energy Transition[END_REF][START_REF]#Envision2030 Goal 7: Affordable and Clean Energy | United Nations Enable[END_REF][START_REF]Paris Agreement[END_REF][START_REF]Clean energy for all Europeans package[END_REF][START_REF]SHARES (Renewables) -Eurostat[END_REF](11)[START_REF]Heating and cooling (facts and figuers)[END_REF][START_REF]District Energy in Sweden[END_REF][START_REF]District Energy in Latvia | Euroheat & Power[END_REF][START_REF]District Energy in the Netherlands[END_REF][START_REF]Energy Policies of IEA Countries -Sweden 2019 Review[END_REF][START_REF]Energy Policies of IEA Countries -France 2016 Review[END_REF][START_REF]Statistics | France -Total Primary Energy Supply (TPES) by source (chart)[END_REF], on the other hand, is better in taking the effect of the earth's internal heat into consideration. However, it is not perfect as it assumes a uniform temperature distribution around the pipe outer surface for its entire length.

𝑅 ′ = 1 2𝜋𝑅 1 ℎ 𝑤 + 𝑙𝑛 | 𝑅 2 𝑅 1 | 2𝜋𝑘 2 + 𝑙𝑛 | 𝑅 3 𝑅 2 | 2𝜋𝑘 3 + 𝑙𝑛 | 𝑅 4 𝑅 3 | 2𝜋𝑘 4 + 1 𝑆𝑘 𝑠 (4-16) 𝑆 = 2𝜋 𝑐𝑜𝑠ℎ -1 ( 𝑍 𝑅 4 ) (4-17) 𝑅 ′ = 1 2𝜋𝑅 1 ℎ 𝑤 + 𝑙𝑛 | 𝑅 2 𝑅 1 | 2𝜋𝑘 2 + 𝑙𝑛 | 𝑅 3 𝑅 2 | 2𝜋𝑘 3 + 𝑙𝑛 | 𝑅 4 𝑅 3 | 2𝜋𝑘 4 (4-18)
where ℎ 𝑤 represents the convective heat transfer coefficient of water, while 𝑘 2 , 𝑘 3 , 𝑘 4 and 𝑘 𝑠 denote the thermal conductivities of the carrier pipe, the insulating material, the outer jacket and the soil, respectively.

The convective coefficient of water is calculated using the Nusselt number (𝑁 𝑢 ) and its thermal conductivity (𝑘 𝑤 ) as shown in equation ( 4 Considering a small fraction of the pipe length 𝑑𝑙, the rate of heat loss (𝑑𝑄) from the water at temperature 𝑇 𝑤 to the ambient at temperature 𝑇 𝑎 can be expressed using equation [START_REF]IEA -International Energy Agency[END_REF][START_REF] Mohammadi | Energy hub: From a model to a concept -A review[END_REF][START_REF]Energy Transition[END_REF][START_REF]#Envision2030 Goal 7: Affordable and Clean Energy | United Nations Enable[END_REF][START_REF]Paris Agreement[END_REF][START_REF]Clean energy for all Europeans package[END_REF][START_REF]SHARES (Renewables) -Eurostat[END_REF](11)[START_REF]Heating and cooling (facts and figuers)[END_REF][START_REF]District Energy in Sweden[END_REF][START_REF]District Energy in Latvia | Euroheat & Power[END_REF][START_REF]District Energy in the Netherlands[END_REF][START_REF]Energy Policies of IEA Countries -Sweden 2019 Review[END_REF][START_REF]Energy Policies of IEA Countries -France 2016 Review[END_REF][START_REF]Statistics | France -Total Primary Energy Supply (TPES) by source (chart)[END_REF][START_REF]Statistics | Sweden -Total Primary Energy Supply (TPES) by source (chart)[END_REF](20)[START_REF] Farhangi | The path of the smart grid[END_REF][START_REF] Del Granado | Synergy of smart grids and hybrid distributed generation on the value of energy storage[END_REF][START_REF] Driesen | Distributed generation: challenges and possible solutions[END_REF].

It should be noted that, depending on where the ambient temperature is taken, either of equations [START_REF]IEA -International Energy Agency[END_REF][START_REF] Mohammadi | Energy hub: From a model to a concept -A review[END_REF][START_REF]Energy Transition[END_REF][START_REF]#Envision2030 Goal 7: Affordable and Clean Energy | United Nations Enable[END_REF][START_REF]Paris Agreement[END_REF][START_REF]Clean energy for all Europeans package[END_REF][START_REF]SHARES (Renewables) -Eurostat[END_REF](11)[START_REF]Heating and cooling (facts and figuers)[END_REF][START_REF]District Energy in Sweden[END_REF][START_REF]District Energy in Latvia | Euroheat & Power[END_REF][START_REF]District Energy in the Netherlands[END_REF][START_REF]Energy Policies of IEA Countries -Sweden 2019 Review[END_REF] or [START_REF]IEA -International Energy Agency[END_REF][START_REF] Mohammadi | Energy hub: From a model to a concept -A review[END_REF][START_REF]Energy Transition[END_REF][START_REF]#Envision2030 Goal 7: Affordable and Clean Energy | United Nations Enable[END_REF][START_REF]Paris Agreement[END_REF][START_REF]Clean energy for all Europeans package[END_REF][START_REF]SHARES (Renewables) -Eurostat[END_REF](11)[START_REF]Heating and cooling (facts and figuers)[END_REF][START_REF]District Energy in Sweden[END_REF][START_REF]District Energy in Latvia | Euroheat & Power[END_REF][START_REF]District Energy in the Netherlands[END_REF][START_REF]Energy Policies of IEA Countries -Sweden 2019 Review[END_REF][START_REF]Energy Policies of IEA Countries -France 2016 Review[END_REF][START_REF]Statistics | France -Total Primary Energy Supply (TPES) by source (chart)[END_REF] is used to calculate the value of the thermal resistance (𝑅 ′ ). This heat loss will in effect reduce the temperature of water in the direction of its flow and it can be rewritten as a function of mass flow rate (𝑚̇) and temperature difference (𝑑𝑇 𝑤 ).

𝑇 𝑤_𝑒𝑛𝑑 = 𝑇 𝑎 + (𝑇 𝑤_𝑠𝑡𝑎𝑟𝑡 -𝑇 𝑎 )𝑒𝑥𝑝 (- 𝐿 𝐶 𝑝 𝑚̇𝑅 ′ ) (4-26)
The term (𝑅′/𝐿) in equation (4-26) corresponds to the total thermal resistance of the pipe which is related to the overall heat transfer coefficient (𝑈) as shown in equation (4-27) [START_REF] Bergman | Fundamentals of heat and mass transfer[END_REF]. Putting (4-27) into gives another form of temperature drop equation as given by (4-28). where 𝑈 𝑅𝑠 is the overall heat transfer coefficient between the water in the pipe and the surrounding cylindrical surface area, 𝐴 𝑅𝑠 , with known temperature profile (either at 𝑅 4 or at Z in Fig. 4-1).

Either of equations or are reported in literature [START_REF] Liu | Modelling, assessment and Sankey diagrams of integrated electricity-heat-gas networks in multi-vector district energy systems[END_REF][START_REF] Liu | Combined analysis of electricity and heat networks[END_REF][START_REF] Shabanpour-Haghighi | An Integrated Steady-State Operation Assessment of Electrical, Natural Gas, and District Heating Networks[END_REF][START_REF] Pirouti | Modelling and analysis of a district heating network[END_REF][START_REF] Fedorov | Parallel implementation of a steady state thermal and hydraulic analysis of heat exchanger[END_REF], but assuming a constant thermal resistance (overall heat transfer coefficient). Such assumptions neglect the variation in the convective coefficient of water as a function of its mass flow rates (which is defined in equations [START_REF]IEA -International Energy Agency[END_REF][START_REF] Mohammadi | Energy hub: From a model to a concept -A review[END_REF][START_REF]Energy Transition[END_REF][START_REF]#Envision2030 Goal 7: Affordable and Clean Energy | United Nations Enable[END_REF][START_REF]Paris Agreement[END_REF][START_REF]Clean energy for all Europeans package[END_REF][START_REF]SHARES (Renewables) -Eurostat[END_REF](11)[START_REF]Heating and cooling (facts and figuers)[END_REF][START_REF]District Energy in Sweden[END_REF][START_REF]District Energy in Latvia | Euroheat & Power[END_REF][START_REF]District Energy in the Netherlands[END_REF][START_REF]Energy Policies of IEA Countries -Sweden 2019 Review[END_REF][START_REF]Energy Policies of IEA Countries -France 2016 Review[END_REF][START_REF]Statistics | France -Total Primary Energy Supply (TPES) by source (chart)[END_REF][START_REF]Statistics | Sweden -Total Primary Energy Supply (TPES) by source (chart)[END_REF]). In this dissertation, the overall heat transfer coefficient is considered as a variable quantity that depends on the mass flow rate. In some case studies, the temperature at the surface of the earth is taken together with equation (4-16) while in the others the soil temperature at the outer most surface of the pipe is considered in combination with equation [START_REF]IEA -International Energy Agency[END_REF][START_REF] Mohammadi | Energy hub: From a model to a concept -A review[END_REF][START_REF]Energy Transition[END_REF][START_REF]#Envision2030 Goal 7: Affordable and Clean Energy | United Nations Enable[END_REF][START_REF]Paris Agreement[END_REF][START_REF]Clean energy for all Europeans package[END_REF][START_REF]SHARES (Renewables) -Eurostat[END_REF](11)[START_REF]Heating and cooling (facts and figuers)[END_REF][START_REF]District Energy in Sweden[END_REF][START_REF]District Energy in Latvia | Euroheat & Power[END_REF][START_REF]District Energy in the Netherlands[END_REF][START_REF]Energy Policies of IEA Countries -Sweden 2019 Review[END_REF][START_REF]Energy Policies of IEA Countries -France 2016 Review[END_REF][START_REF]Statistics | France -Total Primary Energy Supply (TPES) by source (chart)[END_REF].

The last group of equations in the thermal model are those that relate the inlet and outlet temperatures at nodes where pipe flows at different temperatures are mixing. They are By defining the heat equation in pu as shown in equation , the relationship between base mass flow rate and temperature references is derived as shown in equation .

𝑃 ℎ(𝑝𝑢) = 𝑚̇( 𝑝𝑢) ∆𝑇 (𝑝𝑢) (4-37)

𝑚̇( 𝑏𝑎𝑠𝑒) = 𝑃 ℎ(𝑏𝑎𝑠𝑒) 𝐶 𝑝 (𝑇 𝑟𝑒𝑓2 -𝑇 𝑟𝑒𝑓1 ) (4-38)
It is worth noting that the power injection equations for both electricity and heat expressed in actual units may not necessarily be in the same form as those expressed in pu.

In analogous to the electrical transmission line parameters, the base values for pipe parameters in the heating network need to be defined as they are part of the temperature drop equations. Let ∆𝑇 𝑤_𝑠𝑡𝑎𝑟𝑡 and ∆𝑇 𝑤_𝑒𝑛𝑑 be the temperature differences of the inlet and outlet water flows of a pipe relative to the ambient (surface at radius 𝑅 𝑠 ) respectively. The pipe is defined with length 𝐿 and flow rate 𝑚̇. Recalling equation (4-28), ∆𝑇 𝑤_𝑠𝑡𝑎𝑟𝑡 and ∆𝑇 𝑤_𝑒𝑛𝑑 can be related by . Equation (4-40) expresses equation in pu. If the base value of the pipe length is selected to be numerically equal to the base value of the mass flow rate (both in SI units) as shown in equation , then the temperature drop equation in pu can be expressed using equation which has the same form as . Another equation which needs a pu conversion is the pressure drop equation, which can be defined as shown in equation . Then, the base value of the hydraulic head can be derived to be equal to the square of the base mass flow rate as shown in equation .

the PV hubs, and active and reactive electric power generations at the slack hub. For the heating network, on the other hand, the unknown variables are: the values of heat power generation at the slack hub; return temperatures at all temp_supply hubs; supply temperatures at all temp_return hubs; mass flow rates to/from each hub (to be solved from the mass flow rates inside each pipe); and the hydraulic heads at all non-slack hubs (both on the supply and return pipe networks of the heating network).

The flow chart of the algorithm that is used to solve all the unknown variables of integrated electricity and heating networks is depicted in Fig. 4-2. The algorithms for all the constituent submodules are developed using MATLAB ® [113]. The description of each of the submodules follows below.

Network Data: This section is where the data for all the networks is fed by the user. Generations, loads, coupling matrices, network parameters and initial values for the load flow problem are defined here.

Initialization:

The data fed by the user is extracted into matrices so that they can be easily manipulated in MATLAB ® .

Maximum iteration:

The maximum iteration (500 in this work) is the limit in which the solution has to converge with a predefined tolerance. If the solution is not converging within this maximum iteration, it means that either more iterations are required or the known parameters are technically impossible. In any case, the program quits displaying a warning message, the maximum absolute value of the error and the iteration number.

Electricity subroutine:

In this subroutine, the electrical mismatches and the Jacobians that are specific to the electricity network are calculated using the latest values obtained from the Newton-Raphson subroutine.

Hydraulic subroutine: In this section, the pressure resistance coefficient and the overall heat transfer coefficient for each pipe are updated. In addition, the heat power and hydraulic head mismatches are computed using the latest values obtained from the Newton-Raphson subroutine.

Thermal subroutine:

The new values of supply and return temperatures at each hub, together with the new values of pipe flow rates are used to update the inlet and outlet temperatures of each pipe. These end-point temperatures are then used to calculate the outgoing temperature at each node. The differences between these recalculated outgoing temperatures and their corresponding values returned from the Newton-Raphson subroutine give the temperature mismatch values. The values of two types of mismatch equations are calculated here: supply temperature mismatches at temp_return hubs and return temperature mismatches at temp_supply hubs.

Newton-Raphson subroutine:

Here, the unknown variables are numerically solved from the system of nonlinear mismatch equations. Newton-Raphson iterative method is popular for its quadratic convergence rate and wider application area that encompasses meshed and highly coupled topologies (see Section 4.5.2 for the details). At each iteration, the correction values for the unknown variables are calculated by multiplying the mismatch vector with the inverse of the Jacobian matrix. If it is not yet converging, the Jacobian matrix is updated and the iteration continues. However, depending on the problem at hand, Quasi-Newton methods can also be used to approximate the Jacobian matrix. Even more simplified methods (like the forward and backward sweep electricity load flow analysis) can be implemented in the case of radial distribution networks.

Convergence:

Compares the absolute values of all mismatches with the predefined convergence tolerance value (10 -5 pu is used in this dissertation).

Pressure profile calculation on the return pipe network. The mass flow rates on the supply and return pipes of a given branch are equal implying equal pressure drop, but opposite in direction. Hence, the hydraulic heads in the return pipe network can be calculated node by node starting from a reference hub or using a separate Newton-Raphson iterative method. reactive powers, respectively. The terms 𝑃 ℎ , 𝑃 𝑒𝑝 and 𝑃 𝑒𝑞 represent the heat, active and reactive power injections into the network, respectively. 𝑃 ℎ-𝑖𝑛 , 𝑃 𝑓-𝑖𝑛 , 𝑃 𝑒𝑝-𝑖𝑛 and 𝑃 𝑒𝑞-𝑖𝑛 are the heat, fuel, active and reactive power inputs into the coupling system, respectively. 𝐿 ℎ , 𝐿 𝑒𝑝 , and 𝐿 𝑒𝑞 denote the heat, active and reactive power demands. 𝐵𝑜𝑖𝑙𝑒𝑟 𝑃ℎ𝑔 , 𝐻𝑃 𝑃ℎ𝑔 and 𝐶𝐻𝑃 𝑃ℎ𝑔 represent the heat power outputs of the gas boiler, the HP and the CHP, respectively. The terms 𝐻𝑃 𝐿𝑒𝑝 and 𝐻𝑃 𝐿𝑒𝑞 represent the active and reactive electric power consumptions of the HP while 𝐶𝐻𝑃 𝑃𝑒𝑝𝑔 and 𝐶𝐻𝑃 𝑃𝑒𝑞𝑔 designate the active and reactive power generations from the CHP, respectively. Further inspecting the results (by computing the inlet and outlet water temperature for each pipe), it is found that the conservation of energy at node 3 is not obeyed in [START_REF] Shabanpour-Haghighi | An Integrated Steady-State Operation Assessment of Electrical, Natural Gas, and District Heating Networks[END_REF]. At node 3, on the supply pipe network, it is observed that there is a mixing of water coming from node 1 and water coming from the source connected at node 3. The two incoming water temperatures are different. Consequently, the outgoing water temperature at node 3 (i.e. the inlet temperature of water in pipe 3-2) is expected to be different from the water temperature of the source connected to node 3 (i.e. 123.54 °C ), unlike to what it is assumed in [START_REF] Shabanpour-Haghighi | An Integrated Steady-State Operation Assessment of Electrical, Natural Gas, and District Heating Networks[END_REF].

Display results:

The accuracy of the two approaches can be further compared by considering the net heat power injection into the network at node 3. The given heat generation and heat demands at node 3 indicates that there is 9 MW of net heat power injection into the network (see Table 5-3). This should be reflected on the nodal mass flow rate injected into the network and its temperatures on the supply and return sides. In [START_REF] Shabanpour-Haghighi | An Integrated Steady-State Operation Assessment of Electrical, Natural Gas, and District Heating Networks[END_REF], the nodal mass flow rate (calculated from conservation of mass) is 25.08 kg/s while the temperature difference between supply and return sides is 74.01 °C. Computing the heat power injection using these values gives 7.76 MW, which contradicts with the given data of 9 MW. On the other hand, the corresponding nodal mass flow rate and the temperature difference found using the proposed methodology are 28.89 kg/s and 74.51 °C, respectively. The heat power injected into the network calculated by these values gives 9.0 MW which exactly matches with the given data. Hence the proposed model gives more accurate results than the one presented in [START_REF] Shabanpour-Haghighi | An Integrated Steady-State Operation Assessment of Electrical, Natural Gas, and District Heating Networks[END_REF].

Table 5-5 summarizes the network parameters and line flow results found using heat transfer coefficients calculated by applying equation (4-16) on the corresponding standard parameters In the legends, 𝐿 𝑒𝑝 , 𝑃 𝑒𝑝𝑔 and 𝐶𝐻𝑃 𝑃𝑒𝑝𝑔 represent the active power demand, sole generation and production from the CHP plants, respectively. Similarly, 𝐿 𝑒𝑞 , 𝑃 𝑒𝑞𝑔 and 𝐶𝐻𝑃 𝑃𝑒𝑞𝑔 denote the reactive power demand, sole generation and production from the CHPs while 𝐿 ℎ , 𝑃 ℎ𝑔 and 𝐶𝐻𝑃 𝑃ℎ𝑔 stand for the heat power demand, sole production and contribution from the CHPs, respectively. 𝑃 𝑓𝑔 represents the fuel power used by the CHP plants that are installed at different hubs.

- For example, the contributions from all the CHP plants in Case 5-B is 169.465 kW of active, 82.988 kvar of reactive and 192.85 kW of heat powers (see . Adding these powers to the locally produced active, reactive and heat powers and subtracting the corresponding loads (refer Table C-3 and Table 5-7(a)) gives the respective distribution losses (5.938 kW, 13.71 kvar and 140.68 kW respectively). These values are exactly the same as the sum of the corresponding power losses calculated by the model on branch basis for each pipe and transmission lines as indicated in Table 5-8(a). The losses account for 1.36%, 5.07% and 25.57% of the respective demands.

It can be noted that hub 2 in Case 5-B was initially set as temp_supply type (see Table C-3 in Appendix C) assuming it could have an excess of heat power to be injected into the network. On the contrary, hub 3 was set as temp_return type assuming it could have a deficit of heat power. But, the results show that there is a deficit of heat power at hub 2 and excess of heat at hub 3 (see the mass flow injections at each hubs in Table 5-7(a)). Hence, the return temperature at hub 2 is fixed to the maximum possible value and the supply temperature is solved. At hub 3, on the other hand, the supply temperature is fixed to the maximum possible value and the return temperature is solved in the load flow. This shows the model's capability to detect if a given hub is acting as a source or a consumer and configure the load flow equations accordingly.

Cases 5-C and 5-D are extensions of Case 5-B with the addition of hub 5, which has no heat and electricity sole generations, but a CHP that takes 55 kW of fuel and produces both heat and electricity. The outputs of the CHP plant are 18.15 kW of active power, 5.5 kvar of reactive power and 13.75 kW of heat power (see Table 5 † The mass flows through the pipe connecting hubs 𝑖 and 𝑗 as observed on the supply pipe network of the DHN.

The active, reactive and heat power losses for Case 5-C (3.37%, 14.04% and 32.2%, respectively) are higher than the losses in Case 5-D (1.39%, 5.17% and 20.45% respectively). One reason for this is the higher number of transmission lines and pipes used in Case 5-C as compared to Case 5-D. The reactive power loss on a transmission line can be negative when the reactive power supplied by shunt capacitive part of the line becomes greater than what is consumed by the inductive part of the transmission line as it is the case for line 2-4 in Case 5-B (see Table 5-8(a)).

As it can be seen from Table 5- Similarly, mass flow rates from each hub into the network are assumed positive if the hub is acting as a heat source and negative otherwise, as shown in Fig. 5-2. The results show that hubs 2, 4 and 5 are acting as a heat sink and hubs 1 and 3 are acting as a heat source (see Table 567). It can also be observed that the sum of all mass flows from each hub into the DHN (indicated by 𝑚 ̇) are zero for all cases which shows the conservation of mass is holding true.

Although Case 5-C and Case 5-D have the same load profiles for all hubs, there is a very significant difference on their hydraulic head profiles (see Fig. 5-4 and Table 567). The hydraulic head, which is the sum of pressure head and altitude, may not need to be exactly equal to the values indicated in Table 5-7. However, the same relative hydraulic head difference between different hubs should be maintained to guarantee the required pipe mass flows in the network. All hydraulic heads on the supply pipe network or all hydraulic heads on the return pipe networks can be increased/decreased by the same amount as far as they remain in the acceptable hydraulic head ranges. These minimum and maximum acceptable hydraulic head limits depend on the altitude and atmospheric pressures of the site where the pipes are buried. Such limits should be guaranteed to avoid boiling of water inside a pipe (minimum limit) and damage of pipes (maximum limit) [START_REF] Frederiksen | District Heating and Cooling[END_REF].

The highest and the lowest hydraulic heads on the supply network side of the DHN in Case 5-C are 135 m and 103.123 m while in Case 5-D they are 135 m and 1.34 m, respectively (see Table 5-7 and Fig. 5-4). It means that the DHN in Case 5-C needs a lower pumping energy than Case 5-D. This could be associated with the topological differences. Due to more parallel paths in Case 5-C, there are lower mass flow rates in each path with a lower pressure difference between the connected hubs as compared to a single radial path in Case 5-D.

The temperature profiles of the heating network, as it can be seen from Fig. 5-4, shows that all the supply temperatures are less than or equal to the maximum of the values specified at temp_supply hubs (i.e. 90 °C at hub 3, Table C-3). Similarly, all the return temperatures are less than or equal to the maximum of the values specified at temp_return hubs (i.e. 35 °C at hubs 4 and 5, Table C-3) * . It can also be observed that all the temperatures are above the temperature of earth's surface (assumed 10 °C). The inlet and outlet temperatures of each pipe in both the supply and return pipe networks can be referred from Table 5-8.

The results show that for the same dimension of pipe, the drop in temperature and the heat loss (in percentage of the amount transported) increase as the mass flow rate decreases. Accordingly, the largest temperature drop is observed on pipe 3-2 of Case 5-C (both on supply and return network). It also shows that the heat loss across this pipe is higher than what is actually transported. In this case, it may be economical to isolate the pipe (more discussion can be found in Chapter 6). On the other hand, less pipe mass flow implies less friction and, as a result, less pumping energy required to keep the mass flows.

The voltage profiles for all hubs are depicted in Fig. 5-5. It shows that, although Case 5-C and Case 5-D have the same load and generation profiles, except at the slack hub, they experience different voltage profiles due to their topological differences. In all cases, the slack hub and hub 3 (which is PV hub) have a voltage magnitude of 1.0 pu and 1.05 pu, respectively, which is in agreement with their definition in Table C-3. Energy Information database [122]. All the heating demands, including domestic hot water, are assumed to be supplied through the DHN. The DHN pipes are assumed to be buried at 0.5 m depth from the surface. The soil temperature at this depth for the 192 hours is taken from AgriMet online weather database [123]. A uniform temperature profile is assumed around all the DHN pipes. Accordingly, equation [START_REF]IEA -International Energy Agency[END_REF][START_REF] Mohammadi | Energy hub: From a model to a concept -A review[END_REF][START_REF]Energy Transition[END_REF][START_REF]#Envision2030 Goal 7: Affordable and Clean Energy | United Nations Enable[END_REF][START_REF]Paris Agreement[END_REF][START_REF]Clean energy for all Europeans package[END_REF][START_REF]SHARES (Renewables) -Eurostat[END_REF](11)[START_REF]Heating and cooling (facts and figuers)[END_REF][START_REF]District Energy in Sweden[END_REF][START_REF]District Energy in Latvia | Euroheat & Power[END_REF][START_REF]District Energy in the Netherlands[END_REF][START_REF]Energy Policies of IEA Countries -Sweden 2019 Review[END_REF][START_REF]Energy Policies of IEA Countries -France 2016 Review[END_REF][START_REF]Statistics | France -Total Primary Energy Supply (TPES) by source (chart)[END_REF] is used to calculate the thermal resistance of each pipe in the thermo-hydraulic model.

Table 5-10 summarizes the location, type, lowest and peak values of the electricity and heat demands. Solar PV, waste incineration CHP, gas fired HOB and HPs are considered to supply the energy demands. These technologies are dispatched in a merit order based on their operating cost. A 200 kW peak solar PV generation plant is connected at hub 5. The electricity production from the solar PV is forecasted for a year using HOMER ® [START_REF] Homer Energy | HOMER -Hybrid Renewable and Distributed Generation System Design Software[END_REF]. The generation profiles for the selected eight days are then extracted and considered as active electric power production at hub 5. All the available PV generation is assumed to be fed into the electricity distribution network.

The waste incineration CHP plant is connected at hub 3 (also referred to as CHP3). The peak rating of the CHP is 70 kWe and 80 kWth with the corresponding efficiencies of 35% and 40%. It runs at a lagging power factor of 0.86. The CHP is assumed to have two equal units.

Considering annual inspection and maintenance periods, the units are scheduled as follows:

 both units run at full capacity during winter days.  only one unit runs during autumn and spring days.  both units are turned off during summer days.

The gas boiler is installed at hub 4 with 45 kWth capacity at 90% efficiency. It is scheduled to supply part of the local heat demand at hub 4 exceeding 100 kWth.

The three HPs, 20 kWe capacity each with a COP of 4.0 running at 0.9 lagging power factor, are considered. Two scenarios are defined based on their location. In Scenario I, the HPs are installed at hubs 2, 4 and 6 with designation HP2, HP4 and HP6, respectively. They are scheduled to supply the local heat demand at the corresponding hub with a maximum limit of 80 kWth. In Scenario II, on the other hand, HP2, HP4 and HP6 are altogether installed at hub 3, represented as HP3. The heat produced by HP3 at each hour is equal to the heat that could have been produced by HP2, HP4 and HP6 collectively in Scenario I. The scheduling of the CHP at hub 3, the gas boiler at hub 4 and the solar PV at hub 5 is the same in both scenarios. Table 5-11 summarizes the technologies considered in the two scenarios. The duration curves for active power generation, as shown in Fig. 5-12, are quite similar in the two scenarios. The main reason is the voltage level (4.16 kV) which is high relative to the electricity demand, the size of the HPs and the length of distribution networks. Further investigation considering lower (0.4 kV) distribution networks is presented in Chapter 6.

The active electricity generated in 192 hours from the CHP and solar PV is about 10.08 MWh and 6.142 MWh, respectively. The total active power demand for the connected loads at Hubs 2, 4 and 6 is 34.959 MWh. In addition, 5.018 MWh of active power is consumed by the HPs. It shows that the generations from the CHP and the solar PV cover only 32.2% of the total demand. The active power losses in the network is 257.052 kWh in Scenario I and 245.32 kWh in Scenario II, which accounts only 0.643% and 0.614% of the total active electricity demand. The slight difference between the two scenarios is due to the HPs in Scenario I being located far from hub 3 where there is active electricity source. The very low distribution loss also manifests the advantage of high voltage level relative to the installed electricity demands.

Although the network is importing active electricity most of the time in both scenarios, it also exports occasionally as indicated by the circles in Fig. 5-12. During these hours, the sum of the electricity generated from the CHP and the solar PV is higher than the total demand plus the loss at the corresponding hour. The total imported and exported active electricity in Scenario I is 24.416 MWh and 403.976 kWh, respectively. In other words, the energy system in Scenario I is importing a net of 24.012 MWh active electricity in 192 hours. Similarly, the energy system in Scenario II is importing 24.405 MWh and exporting 404.201 kWh of active electricity with a net import of 24.001 MWh. Like the active electricity, the duration curves of reactive power generation for the two scenarios, as shown in Fig. 5 To compare the two scenarios in terms of operating cost, the following price data are considered. The waste incineration CHP plant and the solar PV are assumed to have zero operational cost. Energy prices of 0.2 €/kWh of active electricity, 0.1 €/kWh of heat from the neighbourhood and 0.113 €/kWh of gas for the boiler are assumed. The pumping cost to overcome the pressure drops across the pipes is also taken into account. In addition, a selling price while exporting is assumed to be equal to 60% of the buying price.

A comparison between radial and meshed networks revealed the possibility of higher distribution losses in the latter case. On the other hand, the pseudo-dynamic simulation case study shows that sizing the generation capacities of coupling technologies based on one type of demand profile could lead to undersized or oversized generation for the other network. The results also showed that the distribution losses and the total operating cost vary depending on the location of heat pumps.

The following chapter illustrates how an optimisation routine can be combined with the integrated load flow model to determine optimal placement of coupling technologies and their economic dispatch. either from a CHP or from a gas boiler. In all of the papers, minimization of losses in the electricity network is considered as a main objective function even for DGs that produce both heat and electricity.

Unlike the DGs in the electricity network, the optimal placement of DGs in the heating network is not studied very well. Marguerite et al. [START_REF] Marguerite | Application of a District Heating Network (dhn) Model for an Ex-Ante Evaluation to Support a Multi-Sources Dh[END_REF] used linear programming to study the optimal dispatch of multi-source DHN while Vesterlund et al. [START_REF] Vesterlund | Optimisation of multi-source complex district heating network, a case study[END_REF] used an evolutionary mixed integer linear programing for thermo-economical optimisation of a multi-source DHN. Neither the optimal placement nor the coupling with the electricity network are considered in their study.

A HP can be taken as a DG for the heating network while it is a distributed load in the electricity network. Thus, the HP will be economical if it is installed near to the electricity generation and near to the thermal load at the same time. However, the electricity generation and the thermal load centre are usually far apart. In that case, both of the electricity and the heating networks should be duly considered to find the optimal placement of the HP. A PSObased algorithm is developed to handle such applications by taking the flexibility and modularity advantages of the integrated load flow model discussed in the previous chapters.

PSO is first proposed in [START_REF] Kennedy | Particle swarm optimisation[END_REF], and is a widely used metaheuristic optimisation algorithm which applies the intelligence of stochastic population of particles to arrive at a global best value [START_REF] Bonyadi | Particle Swarm Optimisation for Single Objective Continuous Space Problems: A Review[END_REF]. It tries to find the global best value in analogy to the way a flock of birds scatter and regroup. The whole group is referred to as a population (swarm) and its individual members are called particles. For a problem with 𝑀 decision variables, the position of each particle is defined by a vector with 𝑀-dimensions, as shown in equation . The new direction and speed of the particle (𝑣 𝑖-𝑛𝑒𝑤 ) is calculated using equation (6-2) which is then applied in equation to find the new position of the particle (𝑥 𝑖-𝑛𝑒𝑤 ).

𝑥 𝑖 = (𝑥 𝑖1 , 𝑥 𝑖2 , … , 𝑥 𝑖𝑀 ) (6-1)

𝑣 𝑖-𝑛𝑒𝑤 = 𝜔𝑣 𝑖𝑜 + 𝑐 1 𝑟 1 (𝑥 𝑖-𝑏𝑒𝑠𝑡 -𝑥 𝑖 ) + 𝑐 2 𝑟 2 (𝑥 𝑔-𝑏𝑒𝑠𝑡 -𝑥 𝑖 ) (6-2)

𝑥 𝑖-𝑛𝑒𝑤 = 𝑥 𝑖 + 𝑣 𝑖-𝑛𝑒𝑤 where 𝜔 is the inertia/damping factor; 𝑣 𝑖𝑜 is the current velocity of particle 𝑖; 𝑐 1 is a selfaccelerating (exploration) factor; 𝑐 2 is a global accelerating (exploitation) factor; 𝑥 𝑖 is the current position of particle 𝑖; 𝑥 𝑖-𝑏𝑒𝑠𝑡 is the best position of the particle in the past; 𝑥 𝑔-𝑏𝑒𝑠𝑡 is the global best position achieved so far; 𝑟 1 and 𝑟 2 are random numbers between 0 and 1.

The maximum iteration and number of particles in a PSO algorithm depend on the nature of the problem, the number of decision variables and the parameters of the PSO itself. Ullmann et al. [START_REF] Ullmann | Comparison of PSO variants applied to large scale optimisation problems[END_REF] suggested fixed values of personal and global acceleration factors (𝑐 1 = 𝑐 2 = 1.494) while Aote et al. [START_REF] Aote | Improved Particle Swarm Optimisation Based on Natural Flocking Behavior[END_REF] argued that better results can be found using a step-by-step variation in the following order: 𝑐 1 decreasing from 2.5 to 0.5 and 𝑐 2 increasing from 0.5 to 2.5. The logic behind such variations is to explore locally in the first iterations and then finetune the global best towards the last iterations. Both approaches are tested for the problem at ii. Random population of N particles is generated in the outer PSO each of them containing whole numbers corresponding to the hubs where the HPs and HOB are located. Placement of multiple units at the same hub is possible.

iii. The position of each particle in the outer PSO is updated according to equation (6-3).

iv.

The equations of the integrated load flow model are configured according to the location of the HPs and the HOB.

v. The inner PSO, shown within the broken rectangle in Fig. 6-1, is used to evaluate the fitness of each particle from the outer PSO loop. A random population of M particles is generated in the inner PSO to find out the best operational sizes of the HPs and HOB together with other generations that have fixed locations (i.e. integrated optimal power flow).

vi. The position of each particle in the inner PSO are updated according to equation (6-3).

vii. Integrated load flow is then computed for each particle in the inner PSO loop which is followed by fitness evaluation using equation .

viii. The fitness values are compared with the particle's best fitness and the global best fitness values in the inner PSO. Accordingly, the particle's best position and the global best position are updated.

ix. Steps (vi)-(viii) are computed for each of the M particles in the inner PSO.

x. If the pre-set maximum iteration of the inner PSO (MaxIter2) is not reached, step (ix) is repeated.

xi. The global best fitness value returned from the inner PSO loop is used as the particle's fitness value in the outer PSO loop.

xii. This value is used to update the particle's best and the global best positions in the outer PSO. The results for both scenarios shows that the optimal placement of HPs is not aligned with the active electricity generation (hubs 1, 3 and 5). They are placed rather at some of the hubs where there is heat demand. In scenario I, 54.6% of the 290.5 kW heat produced and 100 % of the 462.38 kW active power generated are transported in their corresponding distribution networks. In the case of scenario II, however, 73.55% of the 291.3 kW heat and 100% of 459.52 kW active power are delivered through the distribution networks. As a result, scenario I has lower heat loss (1.6% of the heat demand) than Scenario II (1.9% of the heat demand).

The active electric power loss, on the other hand, is higher (3.9% of the demand including HPs' consumption) for Scenario I than Scenario II (3.2% of the demand including HPs' consumption).

The pumping power required for circulation in Scenarios I and II are 0.12 kW and 0.16 kW respectively, which are negligible relative to the total electricity demand. It is also found that only 0.05% and 0.07% of this pumping energy is used to overcome the friction and local pressure losses in the network, respectively, while the remaining is used to keep the pressure difference at the consumer substations. The costs of pumping energy for Scenarios I and II are 0.024 €/h and 0.032 €/h, which are less significant relative to the corresponding cost of heat loss (0.46 €/h and 0.54 €/h respectively). This is in agreement with the findings of [START_REF] Vesterlund | Optimisation of multi-source complex district heating network, a case study[END_REF] but in disagreement with the findings of [START_REF] Wang | Optimizing the District Heating Primary Network from the Perspective of Economic-Specific Pressure Loss[END_REF]. The possible reasons could be the difference in the pipe parameters (diameter, length etc.) and operational strategies. The operating cost of the whole network is 68.3 €/h in Scenario I. This value is reduced to 67.68 €/h in Scenario II which shows nearly 1% saving.

Fig. 6-3 depicts the temperature and mass flow profiles at different hubs. A positive mass flow implies that the hub is acting as a source of heat. Hub 6, for example, is acting as a source in Scenario I and as a consumer in Scenario II. When the mass flow from the hub is closer to zero, the water gets more time to lose its heat and, as a results, its temperature gets closer to the soil temperature (for example the return temperature at hub 2 in Scenario I). If the mass flow is zero, as it is the case for hub 5 in both scenarios, the steady state temperature of water on both of the return and supply pipes becomes equal to the soil temperature.

The pipe mass flows on the supply pipe network of the DHN are shown in Fig. 6-4. All the flows are within the acceptable range of 9.3 kg/s defined in Section 6.3.2. Negative flow for Pipe 4-6 implies that the flow is opposite to the positive flow direction initially assumed (which was from hub 4 to hub 6).

The heat loss in the distribution network for Scenarios I and II are 9.35% and 9.69% of the heat demand, while the active electricity losses are 5.39% and 3.17% of the electricity consumption, respectively. Compared with Case 6-A, the losses are higher which is partly due to longer branches and larger percentage of power flows in the networks.

Table 6-4 also shows that, although the increase in the heat loss from scenario I to II is not significant, the reduction in the electricity loss is tremendous, which is 41.2%. This can be explained by considering the way HPs are installed. In Scenario I, the HPs are placed at the farthest nodes relative to the CHP plant. Nodes 8, 23 and 29 belong to the 0.4 kV distribution network while node 25 belongs to the 4.16 kV distribution network. On the other hand, all of the four HPs are installed inside or very close to the 4.16 kV network in Scenario II. The current drown by the HPs at each of the 0.4 kV nodes is equal to the amount of current they could draw at 4.16 kV multiplied by the transformer ratio (10.4 times in this case). The higher the current flowing, the more the losses it causes in the distribution networks for the same transmission line parameters.

The pumping power required for circulation in Scenarios I and II are 3.79 kW and 5.3 kW respectively, which are again negligible relative to the total electricity demand. Unlike to Case 6-A, about 50.73% and 59.06% of this pumping energy are used to overcome the friction and local pressure losses in the network, respectively, both of which are significant. However, the costs of the pumping energy are still insignificant when they are compared to the cost of heat lost (2.74% and 3.71% considering the importing prices of heat and electricity). The system operating cost is 404.5 €/h for Scenario I while it is 384.65 €/h for Scenario II which shows a reduction of about 5%. It means that a relatively larger saving can be achieved for larger networks by applying an integrated optimisation rather than considering the heating network alone.

It is important to note that nodes 9 and 10 are as optimal as node 8 to place a HP due to the symmetry in the network dimensions and the connected loads at these nodes (see Table 6-3 and Fig. 6-6). This applies to other symmetrical nodes too. The more asymmetry in the network, the faster the PSO gets to the global optimum. Hubs with no heat consumption and production, for example 13, 19 and 26, can be treated as either sources or consumer in the numerical computation depending on the direction of mass flow obtained from the load flow solution. The nodal mass flow is always zero at these hubs. However, the temperature can vary depending on the assumption. In any case, the heat power injected into the network is zero. The mass and temperature balances at the mixing nodes are also guaranteed. Hubs 13 and 26, for instance, are acting as heat sources in Scenario I, and their supply temperatures are fixed at 70 °C. Their return temperatures are then determined by the network, and are equal to the outgoing return temperatures of the corresponding nodes.

The results show that placement and sizing of HPs only from the heating network point of view could overload and increase the loss in the electricity distribution network, especially when HPs are installed at low voltage distribution networks. It could also lead to unacceptable voltage profiles. Such issues are avoided using an integrated optimisation that considers the constraints and distribution losses in both networks. The integrated optimisation has also resulted in a lower operating cost for the overall system and a better voltage profile in the electricity distribution network.

The thermo-economic optimisation further revealed that the pumping cost may sometimes counter balance the advantages of lower supply temperatures. Hence, it is worthy to consider the total cost of pumping energy in thermo-economic optimisations. On the other hand, lowering the return temperature of the consumers (on the primary side) is found to be always economical.

In the following chapter, thermo-economic optimisation of coupled electricity and heating networks having both electrical and thermal load profiles will be presented. Time-series variation of load profiles and generation resources will also be taken into account. The overall simulation is conducted on an hourly basis. Hence, the temperature of water leaving the storage is required not to vary during each hour. This constraint is handled by limiting the total hot (cold) water mass leaving the storage in each hour to be less or equal to the water mass in each layer. In addition, after following a natural stratified heat flow during each hour, the water in each layer is fully mixed at the end of each hour. The parameters that are used to model a three layer stratified thermal storage water tank are summarised in Table 7-1. Assuming an underground cylindrical water tank and a uniform soil temperature distribution around its surface, the optimal radius of a given storage volume (𝑉 𝑠𝑡𝑜𝑟𝑎𝑔𝑒 ) that minimizes the surface heat loss can be calculated using equation .

𝑍 1 𝑍 2 𝑍 3 (𝑇 𝑠 , 𝑚̇) (𝑇 𝑟 , 𝑚̇) (𝑇 1 , 𝑚 1 ) (𝑇 2 , 𝑚 2 ) (𝑇 3 , 𝑚 3 ) 𝑍 1 𝑍 2 𝑍 3 (𝑇 𝑠 , 𝑚̇) (𝑇 𝑟 , 𝑚̇) (𝑇 1 , 𝑚 1 ) (𝑇 2 , 𝑚 2 ) (𝑇 3 , 𝑚 3 ) (a
𝑅 𝑠𝑡𝑜𝑟𝑎𝑔𝑒 = √ 𝑉 𝑠𝑡𝑜𝑟𝑎𝑔𝑒 2𝜋 ⁄ 3 (7-2)
The state of charge of the thermal storage can be expressed in terms of the water temperature in the three layers. During the charging process, the flow inside the storage is downward from top to bottom (Fig. 7-2(a)). The 'new' hot water inflows at the top layer (if 𝑇 𝑠 is higher than 𝑇 1 ) while the 'older' cold water leaves from the bottom layer. The water flowing from layer 2 to layer 3, for example, is 'new' for layer 3, but old (from previous time) for layer 2.

It means that the temperature in layer 3 should be calculated before updating the temperature in layer 2.

The equations governing the thermal interaction between the layers and with the external environment are adapted from the TRNSYS ® [148] model. The temperatures at the three layers during the charging phase are updated using equations (7-3) -(7-5). -4)

𝑇 1(𝑡) = 1 𝐶 𝑝 𝑚 1 {𝐶 𝑝 𝑚̇∆𝑡(𝑇 𝑠 -𝑇 1(𝑡-∆𝑡) ) - 𝐴 𝑐 𝑘 𝑤 (𝑇 1(𝑡-∆𝑡) -𝑇 2(𝑡-∆𝑡) )∆𝑡 0.5(𝑍 1 + 𝑍 2 ) -𝐴 𝑠1 𝐾 𝑠 (𝑇 1(𝑡-∆𝑡) -𝑇 𝑎(𝑡) )∆𝑡 + 𝐶 𝑝 𝑚 1 𝑇 1(𝑡-∆𝑡) } (7-5)
For 𝑇 2 < 𝑇 𝑠 ≤ 𝑇 1 , the incoming flow bypasses layer 1 and flows directly into layer 2. In that case, the 𝑚̇ term in equation (7-5) is replaced by zero and the term 𝐶 𝑝 𝑚̇∆𝑡 (𝑇 1(𝑡-∆𝑡) -𝑇 2(𝑡-∆𝑡) ) in the equation (7-4) becomes 𝐶 𝑝 𝑚̇∆𝑡 (𝑇 𝑠 -𝑇 2(𝑡-∆𝑡) ). Similarly if 𝑇 3 < 𝑇 𝑠 ≤ 𝑇 2 , the inflow will directly get into layer 3. Hence, the 𝑚̇ terms in equations (7-4) and (7-5) are replaced with zeros and the term 𝐶 𝑝 𝑚̇∆𝑡 (𝑇 2(𝑡-∆𝑡) -𝑇 3(𝑡-∆𝑡) ) in equation (7-3) is changed into 𝐶 𝑝 𝑚̇∆𝑡 (𝑇 𝑠 -𝑇 3(𝑡-∆𝑡) ).

The temperatures at the three layers of the storage during the discharging phase are updated using equations (7-6) -(7-8). technologies (such as HP, CHP, gas boiler, wind turbine, thermal storage), time-series heat and electricity demands, electricity price signal, forecasts of renewable sources and the capacity constraints of the components in the system. Based on the input data, the integrated load flow problem is formulated and the decision variables are identified. This is followed by the definition of the PSO structures. The main difference between Fig such as voltage, mass flow rate, temperature, pressure, heat loss, electricity loss, slack hub generations and amount of power flow through each pipe and transmission line.

Could the location of coupling technologies affect the performance of the overall energy system?

Multiple-hours load flow studies are also conducted in the form of a pseudo-dynamic simulation on coupled electricity and heating networks consisting of a CHP, solar PV, gas boiler and HPs. It has been observed that the sizing of generation and coupling technologies based on only one type of demand profile could lead to undersized or oversized generation capacity for the other network. Comparative studies also showed that the HPs' location affects the distribution losses and operating parameters of both distribution networks.

How can the optimal placement of the coupling technologies be determined? How can they be dispatched optimally without compromising the constituent distribution networks?

For economical dispatch of generation units and coupling technologies, with the full consideration of the distribution networks' parameters, a metaheuristic optimisation algorithm is developed based on the PSO technique. It is applied to find the optimal location and sizes of distributed generation and coupling technologies. Two management strategies are compared: the first strategy takes only the heating network and neglects the electricity network while the second strategy considers both networks. The results show that optimising the location and sizes of HPs only from the heating network point of view could lower the heat loss. But, it is demonstrated that such management strategy may cause overloading and more losses in the electricity distribution network. It could also lead to unacceptable voltage profiles. The problem gets worse especially when large HPs are installed in a low voltage distribution network. On the other hand, an integrated optimisation that considers the constraints and distribution losses in both networks resulted in a lower operational cost for the overall system and a better voltage profile in the electricity distribution network.

The circulation pumps are used to overcome two types of pressure losses. The first corresponds to the frictional losses in the pipes while the second is related to the minimum pressure difference required at the substations (on primary side). Their relative magnitude depends on the size and operating strategy of the heating network. It is observed that the electricity consumption of circulation pumps is generally not significant compared to the total electricity demand. However, as the total cost of pumping energy can compete with the cost of heat loss in the network, especially if there are very cheap heat sources, the pumping cost should be taken into account in the thermo-economic optimisations.

A step-by-step approach of identifying and isolating lossy branches of a meshed DHN is demonstrated using an integrated optimal power flow. Such approach helps to operate the DHN efficiently during off-peak seasons (hours). The thermo-economic optimisation results further revealed that lowering the supply temperature does not necessarily result in a lower operational cost. The pumping cost may sometimes counter balance the advantages of lower time while the heating sector is still monopolised dominantly by municipalities or few companies. The level of data communication and transparency between the two markets is very limited. While integration of the energy sectors is given much emphasis, the market integration is not given equal attention. For instance, there are supports (e.g. EU directives) to install large HPs in the heating sector. However, this might necessitate upgrading the electricity distribution network from where the HPs get power. As there is no such binding regulation/framework to share the additional investment cost, such projects may get delayed. Absence of common platform is one of the contributing factors.

The developed tool can be used to setup the missing platform for stakeholders (such as DSOs, generation operators, municipalities, investors and policy makers) to come together and investigate the technical, environmental and financial aspects of an integrated energy systems. The platform also opens a market-oriented research area by serving as a decision support tool in the process of designing a suitable structure for the future energy market. The simulation and optimisation results can be used as input in developing business models for investors and in identifying incentives for policy makers. That way, the synergies can be realized faster, and the energy transition towards the carbon-neutral energy system can be accelerated.

Appendix B. Equations for calculating elements of the Jacobian matrix

From the electrical and thermo-hydraulic models, expressed in per-unit, the following equations are derived. The superscript (subscripts) 𝑖 and 𝑗 refers to hub numbers and their corresponding points of interconnections (i.e. nodes/busses) while the subscripts 𝑚 and 𝑛 indicate the number of those pipes in their respective order (based on the layout of the network, they are automatically extracted). Equations (B-1) to (B-8) hold true for the electrical model. Let's define a nodal mass flow rate injected from hub 𝑖 into the network be 𝑚̇𝑖. Then, equations (B-9) to (B-13) hold true for the thermo-hydraulic model. Consider the Jacobina matrix in equation (B-14). The first three rows refer to the coupled mismatch equations while the last three refer to the uncoupled mismatch equations discussed in Section 4.5.1. Accordingly, the elements of the Jacobian matrix can be calculated using equations (B-15) to (B-50).

- (C-1) 𝐶 ℎ𝑢𝑏1 = 1 0 0 0.3 0 1 0 0.1 0 0 1 0.36 ; 𝐶 ℎ𝑢𝑏2 = 1 0 0 0.25 0 1 0 0.08 0 0 1 0.4 ; 𝐶 ℎ𝑢𝑏3 = 1 0 0 0.27 0 1 0 0.15 0 0 1 0.3 ;

  -4 et Fig. 6-5 (petite structure) et Fig. 6-8 et Fig. 6-10 (grande structure).

  Fjärrvärmenät; Elektricitetsdistributionsnät; Utökad strategi för energihub; Integrerad belastningsstudie; Integrerat optimalt effektflöde; Energinät med flera bärare; Multienergisystem; Optimal placering; Termisk lagring. Preface This thesis is based on a research performed within the framework of Erasmus Mundus Joint Doctorate SELECT+ program 'Environomical Pathways for Sustainable Energy Services' coordinated by the KTH Royal Institute of Technology, Stockholm, Sweden. The research was part of the intake project 'Opportunity Analysis of Decentralized Electricity Grid and District Heating Interaction'. It was funded with the support from the Education, Audiovisual and Cultural Executive Agency (EACEA) (FPA-2012-0034) of the European Commission.

Fig. 1 - 1 . 2 Fig. 1 - 2 . 3 Fig. 1 - 3 . 4 Fig. 1 - 4 . 4 Fig. 1 - 5 . 5 Fig. 1 - 6 . 9 Fig. 2 - 1 . 12 Fig. 2 - 2 . 16 Fig. 2 - 3 .Fig. 3 - 1 . 24 Fig. 3 - 2 . 25 Fig. 4 - 1 . 33 Fig. 4 - 2 . 40 Fig. 5 - 1 . 47 Fig. 5 - 2 . 49 Fig. 5 - 3 . 53 Fig. 5 - 4 . 58 Fig. 5 - 5 . 59 Fig. 5 - 6 . 60 Fig. 5 - 7 . 62 Fig. 5 - 8 . 62 Fig. 5 - 9 . 63 Fig. 5 - 63 Fig. 5 - 65 Fig. 5 - 12 . 66 Fig. 5 - 13 . 66 Fig. 5 - 14 . 67 Fig. 5 - 15 . 68 Fig. 5 - 16 . 68 Fig. 5 - 17 . 69 Fig. 6 - 1 . 75 Fig. 6 - 2 . 79 Fig. 6 -Fig. 6 - 4 .Fig. 6 -Fig. 6 - 6 .Fig. 6 -Fig. 6 - 8 .Fig. 6 - 9 .Fig. 6 - 10 .Fig. 6 -Fig. 6 -Fig. 7 - 3 .Fig. 7 - 4 .Fig. 7 - 5 .Fig. 7 - 6 .Fig. 7 - 7 .Fig. 7 - 8 .Fig. 7 - 9 .Fig. 7 -Fig. 7 -Fig. 7 - 12 .Fig. 7 - 13 .Fig. 7 - 14 .

 112123134144155169211222162331243225413342405147524953535458555956605762586259635635655126651366514675156851668517696175627966466666869610667374757677787977712713714 Fig.1-1. Trends in the share of renewables in the EU28 energy mix (data taken from[START_REF]SHARES (Renewables) -Eurostat[END_REF]). ............................... 2 Fig.1-2. Share of renewables in the heating and cooling sector in EU28 by 2018 (data taken from[START_REF]SHARES (Renewables) -Eurostat[END_REF]). .......... 3 Fig.1-3. Primary energy supply mix in France and Sweden by 2016 (data taken from[START_REF]Statistics | France -Total Primary Energy Supply (TPES) by source (chart)[END_REF][START_REF]Statistics | Sweden -Total Primary Energy Supply (TPES) by source (chart)[END_REF]). ........................ 4 Fig. 1-4. Electricity generation mix in terms of fuel sources by 2016: (a) France and (b) Sweden (data taken from [20]). ..................................................................................................................................................... 4 Fig. 1-5. Heat generation mix in terms of fuel sources by 2016: (a) France (b) Sweden (data taken from [20]). . 5 Fig. 1-6.Organization of the thesis. ...................................................................................................................... 9 Fig. 2-1. Categories of MES modelling approaches reported in literature. ..........................................................12 Fig. 2-2. Example of an energy hub where electricity, natural gas and district heat are interacting (adapted from [59]). The coupling devices are (1) transformer, (2) electrical battery, (3) CHP, (4) heat exchanger, (5) thermal storage and (6) absorption chiller. ..........................................................................................16Fig. 2-3. General scheme of a highly coupled MES. ............................................................................................20 Fig. 3-1. Structure of a general and flexible MES simulation and optimisation model. .......................................24 Fig. 3-2. Different parts of an extended energy hub .............................................................................................25 Fig. 4-1. Buried heating network pipes (a) Layout of supply and return pipes, (b) cross-section of each pipe. ...33 Fig. 4-2. Flow chart of integrated load flow analysis for coupled electricity and heating networks. ...................40 Fig. 5-1. Examples of energy hubs and their EEH representations showing the interactions between heat, electricity and fuel energy carriers. (a) a hub with a HOB, (b) a hub with a CHP, (c) a hub with a HP. .......................................................................................................................................................47 Fig. 5-2. Network configuration of case studies 5-A to 5-D ................................................................................49 Fig. 5-3. Hubs' contribution in the generation and consumption of electricity, heat and fuel -Cases 5-B to 5-D. .........................................................................................................................................................53 Fig. 5-4. Temperature and hydraulic head profiles after load flow solution -Cases 5.B to5-D...........................58 Fig. 5-5. Voltage profiles of electricity network for Cases 5-B to 5-D. ...............................................................59 Fig. 5-6. The layout of the heating and electricity network considered for Case 5-E. ..........................................60 Fig. 5-7. Duration curves of (a) the heat demand and (b) fuel consumption -Case 5-E. .....................................62 Fig. 5-8. Duration curves for active power demand: (a) Scenario I and (b) Scenario II -Case 5-E ....................62 Fig. 5-9. Duration curves for reactive electric power demand: (a) Scenario I and (b) Scenario II -Case 5-E. ....63 Fig. 5-10. EEH representations showing the interactions between heat, electricity and fuel energy carriers: (a) a hub with a HOB and a HP, (b) a hub with a CHP and HP. ..................................................................63 Fig. 5-11. Duration curves of heat power generated: (a) Scenario I and (b) Scenario II -Case 5-E. ...................65 Fig. 5-12. Duration curves of active electricity generation: (a) Scenario I and (b) Scenario II -Case 5-E. .........66 Fig. 5-13. Duration curves of reactive electricity generation: (a) Scenario I and (b) Scenario II -Case 5-E.......66 Fig. 5-14. Comparision in terms of pumping cost and total operating cost: (a) Scenario I and (b) Scenario II -Case 5-E. ..............................................................................................................................................67 Fig. 5-15. Temperature and mass flow profiles at hub 1: (a) Scenario I and (b) Scenario II -Case 5-E. ............68 Fig. 5-16. Voltage magnitude and angle at hubs 5 and 6: (a) Scenario I and (b) Scenario II -Case 5-E. ............68 Fig. 5-17. End point water temperatures of pipe 2-3: (a) Scenario I and (b) Scenario II -Case 5-E. ..................69 Fig. 6-1. Flow chart of the nested PSO applied for optimal placement and sizing of HPs and a HOB -Cases 6-A and 6-B. ...............................................................................................................................................75 Fig. 6-2. One-line diagram of the heating and electricity network considered in Case 6-A. ................................79

∑ 1 ( 4

 14 (𝑎𝑙𝑙 𝑚𝑎𝑠𝑠 𝑓𝑙𝑜𝑤𝑠 𝑖𝑛𝑡𝑜 𝑡ℎ𝑒 𝑛𝑜𝑑𝑒) = ∑(𝑎𝑙𝑙 𝑚𝑎𝑠𝑠 𝑓𝑙𝑜𝑤𝑠 𝑙𝑒𝑎𝑣𝑖𝑛𝑔 𝑡ℎ𝑒 𝑛𝑜𝑑𝑒) (4-8) 𝐻 𝑗 -𝐻 𝑖 = (𝑠𝑖𝑔𝑛)𝑘 𝑖𝑗 |𝑚̇𝑗 𝑖 | 𝛽 = 𝑘 𝑖𝑗 𝑚̇𝑗 𝑖 |𝑚̇𝑗 𝑖 | 𝛽-

  𝑑𝑄 = -𝐶 𝑝 𝑚̇ 𝑑𝑇 𝑤 (4-24) Equation (4-25) can be then obtained by equating equations (4-23) and (4-24) and rearranging the terms. Then it is solved for temperatures at the ends of the pipe as shown in equation (4-26).

Fig. 4 - 2 .

 42 Fig.4-2. Flow chart of integrated load flow analysis for coupled electricity and heating networks.

1 .

 1 coupling equation governing an energy hub without any coupling technology is given by equation (Examples of energy hubs and their EEH representations showing the interactions between heat, electricity and fuel energy carriers. (a) a hub with a HOB, (b) a hub with a CHP, (c) a hub with a HP. In the figure, 𝑃 𝑓𝑔 , 𝑃 ℎ𝑔 , 𝑃 𝑒𝑝𝑔 and 𝑃 𝑒𝑞𝑔 terms indicate the local generations of fuel, heat, active and

3 .

 3 Hubs' contribution in the generation and consumption of electricity, heat and fuel -Cases 5-B to 5-D.

- 6 )

 6 which are insufficient to supply the demands (90 kW, 40 kvar and 150 kW respectively) of the hub. The deficits and the additional losses in the connecting branches led to more generations at the slack hubs of Case 5-C and Case 5-D, as compared to Case 5-B.

8 ,

 8 mass flows through pipe 3-2 in Case 5-B, pipe 4-5 in Case 5-C and pipe 3-2 in Case 5-D are negative implying that the directions of flow in these pipes are opposite to what are assumed in Fig 5-2.

Fig. 5 - 4 .

 54 Fig. 5-4. Temperature and hydraulic head profiles after load flow solution -Cases 5.B to5-D.

Fig. 5 -

 5 Fig. 5-7 to Fig. 5-9 show the duration curves of the total consumption of each energy carrier.The numbers in the legends indicate the corresponding hubs. Fig.5-7(a) shows the duration curve of heat demand while Fig.5-7(b) depicts the fuel consumption pattern for the gas boiler and the CHP plant. Both of them are plotted using the total heat demand as a reference. The fuel consumption of the gas boiler follows the heat demand of hub 4 (indicated as load4) exceeding the 100 kW threshold, but the CHP scheduling is not strictly following the hourly total heat demand profile as it is scheduled on seasonal bases. The duration curves shown in Fig.5-8 and Fig.5-9 consist of the active and reactive power consumptions of the coupling technologies, in addition to the residential and commercial electrical loads. It should be noted that the total active and reactive power demands are the same in both scenarios. The only difference is the location of the HPs.

Fig. 5 - 12 .

 512 Fig. 5-12. Duration curves of active electricity generation: (a) Scenario I and (b) Scenario II -Case 5-E.

Fig. 5 -

 5 13. Duration curves of reactive electricity generation: (a) Scenario I and (b) Scenario II -Case 5-E.

  xiii. Steps (iii)-(xii) are computed for each of the N particles in the outer PSO. xiv. If the pre-set maximum iteration of the outer PSO (MaxIter1) is not reached, step (xiii) is repeated. xv. The best location of the HPs and HOB is contained in the global best position of the outer PSO while their economical dispatch is contained in the corresponding global best position of the inner PSO.

Fig. 6 -

 6 Fig. 6-7 shows the nodal mass injection from each hub and the corresponding temperature profiles on the supply and return pipes. The nodal mass flows for both Scenarios lookalike for all nodes except where the HPs are installed (hubs 7, 8, 11, 16, 23, 25 and 29 in Fig. 6-7(a)).

2 .

 2 Three layer stratified model of a thermal storage: (a) charging and (b) discharging.

3 .

 3 . 7-3(a) and Fig. 7-3(b) lies on how the respective PSOs are implemented. Start End Determine the number of decision variables (Nvar) for each hour Yes Input data about network topology, coupling technologies, demand and renewable generation forecasts Extract the decision variables for hour h Run the integrated load flow for hour h Particle i fitness = sum (hourly operating costs); i, ; set h = 1 Randomly initialize all the particles; set iter = 1 Calculate the hourly operating cost considering penalty for any violation of constraints Define the PSO structure: Maximum iteration, N particles each with 24*Nvar decision variables Is h < 24 ? initialize the particles and set iter = 1 Calculate the hourly operating cost considering penalty for any violation of constraints Define the PSO structure: Maximum iteration, N particles each with Nvar decision variables Is i < N ? No iter = iter + 1 Run the integrated load flow Particle i fitness = hourly operating costs; update particle s best position Is h < 24 ? Flow chart of thermo-economic optimisation: (a) with a thermal storage (MPC) and (b) without a thermal storage.

  𝑠(𝑖) -𝑇 𝑟(𝑖) ) 𝑖𝑓 𝑡ℎ𝑒 𝑝𝑖𝑝𝑒 𝑓𝑙𝑜𝑤 𝑖𝑠 𝑙𝑒𝑎𝑣𝑖𝑛𝑔 𝑛𝑜𝑑𝑒 𝑖 -(𝑇 𝑠(𝑖) -𝑇 𝑟(𝑖) ) 𝑖𝑓 𝑡ℎ𝑒 𝑝𝑖𝑝𝑒 𝑓𝑙𝑜𝑤 𝑖𝑠 𝑔𝑜𝑖𝑛𝑔 𝑖𝑛𝑡𝑜 𝑛𝑜𝑑𝑒 𝑖 0 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (B-13)

  modéliser le MES soit décrite dans les chapitres 3 et 4, des méthodologies supplémentaires spécifiques au flux de puissance optimal et au MPC sont incluses dans les chapitres 6 et 7, respectivement. Les principales conclusions sont ensuite résumées au chapitre 8.

4. Quelle est l'importance du stockage thermique dans l'exploitation des synergies ? Comment la gérer de manière optimale en présence d'énergies renouvelables intermittentes et de signaux de prix variables tout en tenant compte de toutes les contraintes techniques ? Cette thèse propose une nouvelle méthodologie de modélisation, simulation et optimisation des réseaux intégrés d'électricité et de chauffage basée sur une approche étendue de hub énergétiques (An Extended Energy Hub -EEH). Des études de cas de démonstration sont considérées pour illustrer la capacité de l'outil à répondre aux questions de recherche cidessus en prenant des paramètres et des conditions aux limites très proches des systèmes réels. Objectifs L'objectif principal de cette thèse est de développer une modélisation, une simulation et une optimisation flexibles qui puissent intégrer les paramètres de fonctionnement des réseaux de distribution dans un système générique multi-énergie (MES), et démontrer son application dans un réseau intégré de chauffage et d'électricité. Il se décline en objectifs spécifiques, axés sur le développement :  d'un modèle général et flexible de flux de charge capable de gérer des MES fortement couplés,  de modèles détaillés de flux de charge intégré régissant les réseaux couplés de distribution d' électricité et de chaleur,  d'un algorithme d'optimisation intégré au modèle de flux de charge MES,  d'un sous -modèle de stockage thermique intégré au modèle MES  d'un algorithme de contrôle prédictif du modèle (MPC) pour une gestion optimale du stockage thermique dans un MES. Les objectifs applicatifs suivants sont également considérés pour démontrer la capacité de l'outil développé et sa nouveauté.  Analyse statique et études des séries temporelles des flux des réseaux couplés d'électricité et de chauffage.  Dispatch économique sur une heure, optimisation thermo-économique et placement optimal des technologies de couplage dans les réseaux intégrés de chauffage et d'électricité.  Optimisation thermo-économique sur des séries temporelles.  Démontrer la capacité de l'outil à exploiter la synergie des réseaux couplés d'électricité et de chauffage en présence de stockage thermique, d'énergies renouvelables intermittentes et d'un signal de prix de l'électricité variable. Structure de la thèse La Fig. 1-6 illustre comment la thèse est organisée en 8 chapitres en lien avec les publications (voir la section « Liste des publications » pour plus de détails) . Le premier chapitre couvre le contexte, l'énoncé du problème et les objectifs de la thèse. Le chapitre 2 est consacré à l'état de l'art de la modélisation d'un MES. Dans ce chapitre, les compléments scientifiques nécessaires sont identifiées. Il est suivi des chapitres 3 et 4 qui détaillent et discutent la méthodologie proposée dans cette thèse. Le chapitre 3 décrit le cadre général appelé « An Extended Energy Hub (EEH) » pour la modélisation des MES. À l'aide de ce cadre, les détails d'un modèle de flux de charge intégré composé de réseaux de chauffage et d'électricité sont examinés au chapitre 4. La méthodologie est ensuite suivie par les chapitres 5, 6 et 7 qui présentent les résultats et les discussions sur les études de cas qui sont utilisées pour illustrer la nouveauté et la capacité de la méthodologie proposée. Des études de cas de flux de charge intégrés sont présentées au chapitre 5 tandis qu'une étude de cas de flux de puissance optimale intégrée est discutée au chapitre 6. Le chapitre 7 décrit comment le stockage thermique peut être intégré dans le modèle et présente l'application du MPC pour exploiter la synergie des réseaux couplés d'électricité et de chauffage. Bien que la principale méthodologie utilisée pour

Tadesse Ayele fonctionnement

  et les contraintes du réseau de distribution ne sont pas étudiés. Ces approches de modélisation peuvent être classées en deux : d'un côté celles qui utilisent des progiciels existants et de l'autre celles qui utilisent leurs propres modèles mathématiques. Les études réalisées à l'aide de progiciels peuvent être regroupées en deux. Le premier sous-groupe se compose des études traitant des détails d'une centrale de co/polygénération et de bâtiments tandis que l'autre sous-groupe couvre la littérature se concentrent uniquement sur le bilan énergétique au niveau du bâtiment, de la ville ou du pays. Concernant les travaux basés sur des développements propres, on trouve des approches modulaires et non modulaires.Les modèles de la deuxième grande catégorie prennent en compte les réseaux de distribution et abordent partiellement ou complètement leurs paramètres physiques et de fonctionnement. Ces modèles sont utilisés pour des études de flux de charge, de dispatch économique et des études de flux de puissance optimales. Un certain nombre de travaux ont proposé un environnement de co-simulation en utilisant les outils existants tandis que d'autres proposent leur propre modèle de réseaux de distribution multi-vecteurs.

2-1). Les travaux de littérature de la première catégorie ne prennent pas en compte le réseau de distribution. Ils se concentrent principalement sur l'évaluation de faisabilité d'un MES composé de différentes technologies énergétiques ou sur l'étude paramétrique détaillée d'une seule centrale électrique/bâtiment. Les pertes du réseau de distribution sont soit négligées, soit représentées en pourcentage fixe. Les paramètres de x | Doctoral thesis by Getnet

Les progiciels de modélisation MES existants utilisent l'équilibre énergétique entre la demande agrégée et la source. Les réseaux de distribution ne sont pas correctement modélisés dans de tels outils. Il existe en revanche des outils de simulation et d'optimisation bien développés, capables de gérer un seul vecteur d'énergie. Très peu de travaux ont tenté de les combiner dans un environnement de co-simulation afin d'intégrer les paramètres de fonctionnement des réseaux de distribution. Il existe également travaux proposant des modèles et des algorithmes avec et sans prise en compte des paramètres de fonctionnement du réseau de distribution. La plupart d'entre eux sont spécifiques à chaque cas, suivant une approche modulaire ou non. Les modèles modulaires ont appliqué le concept de hub énergétique sur la base des flux de puissance du réseau. Cependant, l'hypothèse de flux unidirectionnels et de modèles hydrauliques isothermes limite leur champ d'application uniquement à un MES composé des seuls réseaux d'électricité et de gaz, sans aucun prosommateur (prosumer).

Le concept de hub énergétique doit être étendu (redéfini) pour ajouter une flexibilité de flux d'énergie bidirectionnels et une capacité à gérer les réseaux de distribution de chaleur. Par ailleurs, les modèles thermo-hydrauliques doivent être suffisamment généraux pour inclure des topologies maillées et suffisamment flexibles pour que tout changement dans la topologie ou la configuration du réseau de chauffage puisse être automatiquement détecté par l'algorithme. L'évolutivité et la modularité de la méthodologie pour intégrer tout vecteur d'énergie supplémentaire ou toute technologie de couplage additionnelle sont également cruciales.

  Les profils de demande, les matrices de génération et de couplage locaux décrivant les technologies énergétiques installées dans chaque hub sont également intégrés dans le modèle par l'utilisateur. Des informations supplémentaires sur les valeurs de référence, les valeurs de base, les constantes, le coût des générations, etc. sont également fournies à ce niveau.A partir de ces données d'entrée, la partie core process permet l'extraction automatique des informations requises pour formuler le problème de flux de charge. Les paramètres décrivant les différents vecteurs d'énergie sont convertis à l'aide d'un système unitaire commun (p.u.). Une distinction est alors faite entre les variables de référence (qui ont des valeurs fixes) et celles qui doivent être résolues. Ensuite, un système d'équations, également appelé problème de flux de charge intégré, est formulé et résolu pour ces variables inconnues. Si l'objectif opérationnel porte sur l'optimisation, le problème de flux de charge sera utilisé comme une contrainte d'égalité parmi d'autres contraintes d'inégalité. Comme le problème de En remplaçant l'équation (3-2) par l'équation (3-1), la version étendue de l'équation de couplage pour un réseau général de MCES peut être formulée comme indiqué dans l'équation. Comme les injections de puissance dans le réseau (qui font désormais partie de l'équation de couplage) sont des fonctions des paramètres du réseau à tous les hubs, elles décrivent le couplage entre les demandes et les productions au niveau de chaque hub.Le point d'interconnexion est l'endroit où l'injection nette d'énergie (c'est-à-dire 𝑃 𝑎 , 𝑃 𝑏 , …, 𝑃 𝑛 sur la Fig.3-2 du hub dans le réseau a lieu pour différents vecteurs d'énergie. Il correspond au terme bus dans le réseau électrique et noeud dans le réseau de chauffage urbain.

	que les variables inconnues sont résolues) si et seulement si les amplitudes de tous les
	décalages sont inférieures à une tolérance acceptable.
	Les réseaux couplés d'électricité et de chauffage sont représentés à l'aide de trois sous-
	modèles : les modèles électrique, hydraulique et thermique. Bien que la représentation
	rectangulaire de paramètres complexes, tels que la tension et l'impédance, puisse être utilisée
	pour la modélisation de réseaux électriques, la représentation hybride (profil de tension
	polaire et impédance rectangulaire) est la plus largement utilisée pour formuler les équations
	La d'équilibre de puissance [99-101]. Les modèles hydrauliques basés sur des équations d'écoulements en conduites aboutissent à un modèle simple avec un taux de convergence
	distance entre le système de couplage et ce point d'interconnexion est supposée très courte aussi rapide (voire meilleur) que les équations nodales et en boucle [94]. Par conséquent,
	par rapport à la longueur des branches reliant deux hubs différents. Cela permet de négliger les éventuelles chutes de tension, de pression et de température entre le système de couplage et le point d'interconnexion. Cependant, cela ne signifie pas qu'ils sont une seule et même chose. Par exemple, si des flux d'eau ayant des températures différentes se mélangent cette méthode est utilisée pour représenter la partie hydraulique du réseau de chauffage. Le modèle thermique est représenté à l'aide de trois groupes d'équations qui régissent l'injection d'énergie thermique, la chute de température le long d'une conduite du fait des pertes en ligne et la conservation de l'énergie lors du mélange de l'eau à chaque noeud.
	à un point d'interconnexion donné connecté à un hub de source de chaleur, les températures Une fois que toutes les équations de couplage à chaque hub et les équations bilans pour de l'eau sortante côté alimentation seront différentes de la température d'alimentation du chaque réseau de distribution sont formulées dans un système normalisé, le système concentrateur hub. Il est donc important de traiter séparément le hub énergétique et le point d'équations (également appelé problème de flux de charge) est résolu comme un problème d'interconnexion afin de modéliser correctement les réseaux thermiques. unique en utilisant la méthode itérative de Newton-Raphson. Les variables d'intérêt dans la La partie réseau d'énergie, représente le reste du système global -le réseau de distribution solution de flux de charge sont l'amplitude de la tension, l'angle de tension, la température
	multi-vecteur (voir Fig. 3-2). Si, pour un vecteur, il y a un excès de puissance à la quantité d'alimentation, la température de retour, la pression hydraulique et les débits massiques de
	requise dans le hub, il sera injecté sur le réseau de ce vecteur et transporté vers les autres hubs chaque hub. À un hub donné, certaines des variables peuvent être connues tandis que les
	du réseau où il y a un déficit. Si le hub présente plutôt un déficit d'un des vecteurs, il autres sont inconnues. L'étude de flux de charge a pour objectif de résoudre les variables
	l'importera du réseau correspondant. Cependant, il convient de noter que le type et la quantité de vecteur d'énergie à fournir/importer vers/depuis le réseau dépendent non seulement des demandes et des générations locales, mais également des types de dispositifs de couplage inconnues à l'aide des variables connues. Selon les variables connues, les hubs peuvent être classés comme indiqué dans le Table 4-1. L'organigramme de l'algorithme développé pour résoudre toutes les variables inconnues du réseau d'électricité et de chauffage est représenté
	disponibles sur un hub donné. En d'autres termes, selon les types de dispositifs de couplage sur la Fig. 4-2. Tout le code est développé à l'aide de MATLAB [113].
	disponibles dans le système de couplage, un type donné de demande locale peut être satisfait
	par une combinaison de puissance entrant dans le concentrateur à partir de différents vecteurs
	d'énergie. S'il y a une surproduction d'un vecteur d'énergie donné à l'intérieur du hub (en
	raison de la conversion entre différents porteurs d'énergie), il sera réinjecté dans le réseau.
	Dans l'étude des flux de charge, seuls certains des paramètres de fonctionnement (variables)
	sont connus et les autres doivent être résolus. Tout calcul de l'injection de puissance basé sur
	des valeurs inexactes des variables inconnues s'écartera des injections de puissance exactes à chaque hub, calculées à l'aide de l'équation (3-3). La différence entre les deux est donnée par l'équation de discordance (mismatch equation). Les équations pour les injections de	et
	l'efficacité de ces dispositifs de couplage. L'équation (3-1) montre la représentation matricielle générale de la relation entre les demandes (m types), les entrées de puissance nette puissance à chaque point d'interconnexion (en fonction des paramètres de fonctionnement) sont combinées avec les équations de couplage (données par l'équation (3-3) pour définir
	dans le système de couplage (n types) et le taux de transfert d'énergie dans les unités de toutes les équations de discordance à résoudre pour déterminer les variables inconnues. Cet
	stockage (n types) à un hub [59]. ensemble d'équations de discordance est appelé le problème de flux de charge d'un MES
	donné. Des méthodes numériques, telles que la méthode de Newton-Raphson, sont utilisées
	pour réduire les amplitudes de non-concordance en modifiant (approximant) de façon
	itérative les valeurs des variables inconnues. La solution est dite convergente (c'est-à-dire

3-1. Il comporte trois parties : input data, core process et outputs. La partie input data correspond aux informations fournies par l'utilisateur parmi lesquelles figurent les paramètres du MES. Les données concernant les réseaux de distribution incluent leurs schémas, les paramètres physiques et les contraintes. flux de charge, la routine d'optimisation est également construite pour se configurer automatiquement avec très peu voire pas de changements de la partie core process, si l'utilisateur modifie les paramètres d'entrée, tels que la topologie du réseau et/ou l'emplacement et les types de technologies énergétiques impliquées. Une fois que l'optimisation ou le problème de flux de charge est résolu, les valeurs unitaires des variables sont reconverties à leurs valeurs réelles. Enfin, soit les solutions optimales ou des résultats de simulation flux de charge en fonction du mérite sont retournés à l'utilisateur (partie outputs de la méthodologie).

Dans la partie core process de la structure proposée à la Fig.

3-1

, une nouvelle approche de modélisation du hub énergétique étendu (EEH) est utilisée pour définir le cadre général qui peut être appliqué à n'importe quel problème de flux de charge intégré régissant tout type de MES. La Fig.

3

-2 illustre les trois parties composant un Hub Energétique dans sa définition étendue proposée dans ce travail: le hub énergétique, les points d'interconnexion et les réseaux d'énergie.

Le hub énergétique comprend différents types de puissances et de demandes locales générées localement ainsi qu'un système de couplage. Le système de couplage est associé à différents types de dispositifs de conversion d'énergie tels que les éoliennes, les panneaux photovoltaïques, les capteurs solaires, les échangeurs de chaleur, les transformateurs de puissance, les chaudières électriques, les chaudières à gaz, les radiateurs électriques, les centrales de cogénération, les PAC, les stockages de chaleur, les stockages de carburant, stockages d'électricité. Les matrices de couplage sont extraites en formulant la relation algébrique entre les demandes, les entrées de puissance nette dans le système de couplage En se référant à la Fig.

3-2

, on peut voir que chacune des puissances d' entrée nettes dans le système de couplage, (𝑃 𝛾_𝑖𝑛 ) termes de l'équation (3-1), est égale à la génération locale correspondante (𝑃 𝛾𝑔 ) moins l'injection dans le réseau (𝑃 𝛾 ) du vecteur d'énergie correspondant (voir équation (3-2)).

Table 5

 5

	-1 to Table 5-4. Les résultats montrent que le modèle proposé
	dans cette thèse est plus précis.
	La contribution de chaque hub et pour chaque vecteur en terme de production et de
	consommation globales des trois	structures de réseaux considérées, est illustrée

  am so grateful to the SELECT+ consortium for giving me a chance to pursue my PhD. I would also like to thank the European Commission, Audiovisual and Cultural Executive Agency (EACEA) for partly funding my study. I truly appreciate all the staffs in the DSEE at IMT Atlantique and in the EKV/EGI at KTH who made my stay at the two institutes comfortable and fruitful. Special thanks to Chamindie Senaratne, for facilitating all the SELECT+ administrative issues and for all the support during my stay at KTH. I would like to express my sincere gratitude to my principal supervisors Prof. Bruno Lacarrière at IMT Atlantique and Prof. Björn Laumert at KTH, and also to my co-supervisor Assoc. Prof. Pierrick Haurant at IMT Atlantique. I am deeply grateful for their immense support, persistent guidance, encouragement and for the trust they put on me throughout the course of this PhD. Moritz Wegener, Mahrokh Samavati, José Fiacro and all of those who are affiliated to the SELECT+ PhD programme who shared me their experiences. I owe my deepest gratitude to Sylvestre Yaovi and Getenet Tesega who encouraged me throughout the journey of this PhD. My special thanks go to Daniel Minilu and Tigist Fetene for their kind hospitability. This dissertation is basically based on the contents of the four journal papers, ([Paper 1], [Paper 2], [Paper 3] and [Paper 4]). Though not directly referred, the extended and revised versions of [Paper 5] and [Paper 6] are also used to complement the thesis. The contribution of the authors in the six papers is presented as follows.[Paper 1] presents the general framework of MES load flow model based on an extended energy hub (EEH) approach. It also discusses the details of an integrated load flow model governing coupled heating and electricity distribution networks. The proposed model is tested and analytically validated by considering radial and meshed topologies of electricity and heating network topologies coupled using multiple combined heat and power (CHP) plants. Also, a comparison is made with a case study from literature. In this paper, the author of this dissertation contributed by designing the conceptual framework of the research, going through the literature review, developing the models, running the simulations, interpreting the results, drawing the conclusions, writing the original draft manuscript, and in the revision and editing process of the manuscript. The research was conducted under the supervision of Prof. Bruno Lacarrière, Assoc. Prof. Pierrick Haurant and Prof. Björn Laumert who participated in the discussions, and reviewed the draft manuscript. Paper 5], the author of this thesis contributed by designing the conceptual framework of the research, going through the literature reviews, developing the models, running the simulations, interpreting the results, drawing the conclusions, writing the original draft manuscripts, and in the revision and editing process of the manuscripts. Prof. Bruno Lacarrière, Assoc. Prof. Pierrick Haurant, Prof. Björn Laumert, and Assoc. Prof. Mohamed Tahar Mabrouk participated in the discussions, and reviewed the draft manuscripts. minimizing the operational cost of coupled heating and electricity networks. It is an invited paper based on its previous version, which was presented at ECOS 2019 international conference, [Paper 6]. In both versions, the author of this dissertation contributed by designing the conceptual framework of the research, going through the literature reviews, developing the models, running the simulations, interpreting the results, drawing the conclusions, writing the original draft manuscripts, and in the revision and editing process of the manuscripts. The work was done under the guidance of Prof Massimo Santarelli who, together with Prof. Bruno Lacarrière, Assoc. Prof. Pierrick Haurant, Assoc. Prof. Mohamed Tahar Mabrouk and Prof. Björn Laumert, participated in the discussions, and reviewed the draft manuscripts.

	Nomenclature | xxxi
	Nantes, France (Dec. 2020) Getnet Tadesse Ayele European Union member countries (including United Kingdom) Author's contributions Nomenclature Abbreviations and acronyms AC Alternating current CHP Combined heat and power COP Coefficient of performance of a heat pump DC Direct current DG Distributed generation DHN District heating network DSO Distribution system operator EEH Extended energy hub EU European Union EU28 HP Heat pump MES Multi-energy system MPC Model predictive control PAC Pompes à chaleur PQ hub A hub where active and reactive power generations are known PSO Particle swarm optimisation PV hub A hub where active power and voltage magnitude are known pu Per-unit Solar PV Solar photovoltaic temp_return hub A hub that imports heat from the network, its return temperature is known temp_supply hub A hub that injects heat into the network, its supply temperature is known TPES Total primary energy supply Latin symbols 𝑎 Transformer tap setting 𝐴 𝑅𝑠 Cylindrical surface area around a pipe at a radius 𝑅 (m 2 ) B Susceptance of transmission line (S) 𝐵𝑜𝑖𝑙𝑒𝑟 𝑃ℎ𝑔 Heat produced by a boiler (W) 𝑪 Coupling matrix 𝐶 𝛿𝛾 An element in a coupling matrix showing how much of the input power of type 𝛾 is converted into the local demand type of 𝛿 𝐶𝐻𝑃 𝑃𝑒𝑝𝑔 Active power produced by a combined heat and power plant (W) 𝐶𝐻𝑃 𝑃𝑒𝑞𝑔 Reactive power produced by a combined heat and power plant (var) 𝐶𝐻𝑃 𝑃ℎ𝑔 Heat produced by a combined heat and power plant (W) 𝐶 𝑙𝑜𝑐 A fraction to take the local pressure losses due to valves and junctions into account 𝐶 𝑝 Specific heat capacity of water (J.kg -1 .K -1 ) Cost of heat locally generated (€/h) Cost of imported/exported heat through the slack hub (€/h) Personal acceleration factor in the PSO algorithm Global acceleration factor in the PSO algorithm Internal diameter of a carrier pipe (m) Infinitesimal length of a pipe (m) Infinitesimal heat loss across an infinitesimal length of a pipe (W) Infinitesimal change of water temperature (K) Internal surface roughness of a pipe (m) A rate of change of stored thermal energy in a hot water storage tank (W) Power from a storage (rate of change of stored energy) (W) friction factor Objective function of an optimization Gravitational acceleration (m.s -2 ) A multivariable system of equations Hydraulic head (m) A vector of unknown hydraulic heads (pu) Base value for hydraulic head (m) Active power consumed by a heat pump (W) Reactive power consumed by a heat pump (var) Heat produced by a heat pump (W) Convective heat transfer coefficient of water (W.m -2 .K -1 ) Current flow (A) Base value for current flow (A) Ampacity of a transmission line (A) Jacobian matrix 𝐶 𝐶 𝑃ℎ𝑔 𝐶 𝑃ℎ(1) 𝑐 1 𝑐 2 𝐷 𝑑𝑙 𝑑𝑄 𝑑𝑇 𝑤 e 𝐸 ̇ℎ 𝐸 ̇𝛾 𝑓 ℱ 𝑔 𝑭 𝐻 𝑯 𝐻 𝑏𝑎𝑠𝑒 𝐻𝑃 𝐿𝑒𝑝 𝐻𝑃 𝐿𝑒𝑞 𝐻𝑃 𝑃ℎ𝑔 ℎ 𝑤 𝐼 𝐼 𝑏𝑎𝑠𝑒 𝐼 𝑚𝑎𝑥 𝑱 𝑘

I want to express my greatest appreciation to Assoc. Prof. Francesco Fuso Nerini, from EGI/KTH, for evaluating the final seminar and giving valuable feedback on the draft thesis that helped me to improve its quality. I am also thankful to Assoc. Prof. Anders Malmquist from EGI/KTH and Fabian Levihn from Stockholm Exergi for their constructive feedbacks on the progresses of this research.

I am grateful to all the jury members for their time in reviewing this thesis and also for their participation in the PhD defence. Many thanks to Prof. Vittorio Verda (Politecnico di Torino), Prof. Bertrand Raison (Université Grenoble Alpes), Prof. Louise Ödlund (Linköping University), Prof. Filip Johnsson (Chalmers University of Technology), Assoc. Prof. Marianne Boix (Université de Toulouse INP-ENSIACET) and Assoc. Prof. Mikael Amelin (KTH Royal Institute of Technology).

I take this chance to thank Assoc. Prof. Mohamed Tahar Mabrouk (DSEE/IMT Atlantique) for the discussions and inputs as a co-author in most of the publications made out of this study. Also to Prof Massimo Santarelli (Politecnico di Torino), I am grateful for his feedback and contribution as a coauthor in two of the publications. I am also thankful to Monica Arnaudo and Johan Dalgren from EGI/KTH for the fruitful discussions we had. I want to show my appreciation to my colleagues and friends who contributed to my study and personal life in one way or the other. Thanks to Imtisal-e-Noor, Manuel Betancourt, Mexitli Eva, The methodology proposed in [Paper 1] is further expanded in [Paper 5] to take hourly variations of demand and generations into account. Pseudo-dynamic simulations on integrated heating and electricity distribution networks consisting of CHPs, heat pumps (HPs) and heat-only-boiler (HOB) are investigated.

[Paper 2] demonstrates how the EEH approach can be integrated with optimisation techniques, such as Particle Swarm Optimisation (PSO), in order to conduct an integrated optimal power flow in coupled heating and electricity distribution networks. It shows how the model can overcome the challenge of having "automatically configurable" model that is crucial in finding the optimal placement of coupling technologies. [Paper 4] further demonstrates the capacity of the proposed modelling and optimisation tool in exploiting the synergy, quantifying the flexibility and capturing key operating parameters of coupled heating and electricity distribution networks in the presence of intermittent renewables and variable price signals. It incorporates a detailed stratified thermal storage model and a PSO-based model predictive control algorithm. In [Paper 2], [Paper 4] and [[Paper 3] introduces exergy and energy analysis as a means of identifying and isolating lossy branches in an electrified district heating network. It also highlights the application of thermo-economic optimisation in 𝑃𝑒𝑙𝑔 Cost of electricity locally generated (€/h) 𝐶 𝑃𝑒𝑙(1)

Cost of imported/exported electricity through the slack hub (€/h) 𝐶 𝑃𝑓𝑔 Cost of fuel used in the system (€/h) 𝑖𝑗 Pressure resistance coefficient for the pipe flow connecting nodes 𝑖 and 𝑗 (m.s 2 .kg -2 ) 𝑘 2

Table 2 -

 2 1: Summary of MES modelling approaches. ........................................................................................21

Table 4 -

 4 1: Types of hubs and the associated known and unknown variables. ...................................................38

Table 5 -

 5 1: Comparison of hub parameters for electricity network -Case 5-A. ................................................50

Table 5 -

 5 2: Comparison of line flows in the electricity network -Case 5-A. .....................................................51

Table 5 -

 5 3 : Comparison of hub parameters for heating network -Case 5-A. ...................................................51

Table 5 -

 5 4: Comparison of branch parameters in the heating network -Case 5-A. ...........................................51

Table 5 -

 5 5: Hub parameters and line flows in the heating network using standard pipe parameters with standard insulation layer -Case 5-A. ...........................................................................................................52

Table 5 -

 5 6: Input and outputs of the CHP plants at different hubs -Cases 5-B to 5-D. .....................................54

Table 5 -

 5 7: Hub-level operating parameters of the electricity and heating networks. ........................................55

Table 5 -

 5 8: Selected branch-level operating parameters in the electricity and heating networks. ......................56

Table 5 -

 5 9: Pipe and transmission line parameters for Case 5-E. .......................................................................60

Table 5 -

 5 10: Summary of the demand profile considered in Case 5-E. ..............................................................61

Table 5 -

 5 11: Summary of the technologies considered in the two scenarios -Case 5-E. ...................................62

Table 6 -

 6

1: Load distribution for Case 6-A. ........................................................................................................79 Table 6-2: Optimal location and dispatches of the HPs and the HOB for Case 6-A. ........................................80 Table 6-3: Load distribution for Case 6-B. ........................................................................................................83 Table 6-4: Optimal location and dispatches of the HP and the HOB for Case 6-B. ..........................................83 Table 6-5: Heat demand at different hubs -Case 6-C. ......................................................................................88

Table 6 -

 6

6: Hub-level results of the optimal power flow (base case) -Case 6-C. ..............................................90

Table 6 -

 6

7: Branch-level results of the optimal power flow (base case) -Case 6-C. .........................................91

Table 6 -

 6

8: Hub-level results of the optimal power flow (reduced topology) -Case 6-C. .................................92

Table 6 -

 6

9: Branch-level results of the optimal power flow (reduced topology) -Case 6-C .............................93

Table 6 -

 6

10: Hub-level results after the optimal placement and sizing of the HPs -Case 6-C. .........................93 Table 6-11: Branch-level results after optimal placement and sizing of the HPs -Case 6-C............................93

Table 6 -

 6

12: Hub-level results of the thermo-economic optimisation -Case 6-C. .............................................94

Table 6 -

 6

13: Branch-level results of the thermo-economic optimisation -Case 6-C. ........................................94

Table 7 -

 7 1: Parameters describing a stratified thermal storage model ..............................................................100

Table 7 -

 7 2: Summary of the optimisation results for the two scenarios in Case 7-A ........................................113 Table C-1: Transmission line and transformer parameters for the electricity networks -Cases 5-B to 5-D ...141 Table C-2: Carrier pipe, insulation and soil parameters for the heating networks -Cases 5-B to 5-D ...........141 Table C-3: Generation and demands of electricity, heat and fuel power at different hubs -Cases 5-B to 5-D ..................................................................................................................................................141Table D-1: Pipe and transmission line parameters for Case 6-C. ....................................................................143 Table D-2: Hub-level results of the optimal power flow after isolating the 1 st lossiest branch -Case 6-C. ....143 Table D-3: Branch-level results of the optimal power flow after isolating the 1 st lossiest branch -Case 6-

C. ..................................................................................................................................................143 Table D-4: Hub-level results of the optimal power flow after isolating the 2 nd lossiest branch -Case 6-C. ...143 Table D-5: Branch-level results of the optimal power flow after isolating the 2 nd lossy branch -Case 6-C ...143

  .......................................................................................................................................................technologies are then discussed in Chapter 6. Chapter 7 illustrates how the tool can be further applied in exploiting the synergy from coupled electricity and heating networks in the presence of intermittent renewables and variable price signals. In addition to the main methodology, which is covered in Chapters 3 and 4, specific methodologies on the optimal placement of coupling technologies and application of MPC are included in Chapters 6 and 7, respectively. The main findings, recommendations for improvement and research outlooks are then summarised in Chapter 8. State of the Art in Modelling Multi-Energy Systems | 19

					Case studies, results
					and discussions
					5. Demonstrations and
					theoretical validation
			Methodology	[Paper 1]
	1. Introduction	2. State of the art in modelling MES	3. General framework to model MES load [Paper 1] flow problem	[Paper 1] networks 4. Integrated load flow model for coupled electricity and heating	7. Applications in the 6. Optimal placement of coupling technologies [Paper 3] [Paper 2]	8. Closure
					presence of intermittent
					renewables and variable
					price signals
					[Paper 4]
			Fig. 1-6. Organization of the thesis.

2.

  -12)[START_REF] Larock | Hydraulics of Pipeline Systems[END_REF]. The Reynolds number, as a function of the water flow rate (𝑚̇) and its dynamic viscosity (𝜇), is defined as shown in equation[START_REF]IEA -International Energy Agency[END_REF][START_REF] Mohammadi | Energy hub: From a model to a concept -A review[END_REF][START_REF]Energy Transition[END_REF][START_REF]#Envision2030 Goal 7: Affordable and Clean Energy | United Nations Enable[END_REF][START_REF]Paris Agreement[END_REF][START_REF]Clean energy for all Europeans package[END_REF][START_REF]SHARES (Renewables) -Eurostat[END_REF](11)[START_REF]Heating and cooling (facts and figuers)[END_REF][START_REF]District Energy in Sweden[END_REF].Due to the transcendental nature of the Colebrook equation (4-12), it is not convenient to use it directly in computer programs. Winning and Coole[START_REF] Winning | Explicit Friction Factor Accuracy and Computational Efficiency for Turbulent Flow in Pipes[END_REF] made a comparative study on the accuracy and computational efficiency of different formulae used to approximate the Colebrook equation and concluded that the Buzzelli equation given by equation[START_REF]IEA -International Energy Agency[END_REF][START_REF] Mohammadi | Energy hub: From a model to a concept -A review[END_REF][START_REF]Energy Transition[END_REF][START_REF]#Envision2030 Goal 7: Affordable and Clean Energy | United Nations Enable[END_REF][START_REF]Paris Agreement[END_REF][START_REF]Clean energy for all Europeans package[END_REF][START_REF]SHARES (Renewables) -Eurostat[END_REF](11)[START_REF]Heating and cooling (facts and figuers)[END_REF][START_REF]District Energy in Sweden[END_REF][START_REF]District Energy in Latvia | Euroheat & Power[END_REF] is the most accurate one. Hence, this equation is used to determine the friction factor for the turbulent and transition intervals (i.e. for 𝑅 𝑒 >2300).

	𝑓 =	64 𝑅 𝑒				(𝑓𝑜𝑟 0 < 𝑅 𝑒 < 2300)	(4-11)
	1 √𝑓	= -2𝑙𝑜𝑔 (	𝜀 3.7𝐷	+	2.51 𝑅 𝑒 √𝑓	)	(𝑓𝑜𝑟 𝑅 𝑒 > 4000)	(4-12)
	𝑅 𝑒 =	4𝑚π	𝐷𝜇				(4-13)
	𝑓 = [𝐵 1 -	𝐵 1 + 2𝑙𝑜𝑔 1 + 𝐵 2 2.18 𝐵 2 𝑅 𝑒	]	-2	(𝑓𝑜𝑟 𝑅 𝑒 > 2300)	(4-14)
	where 𝐵 1 =	0.774 ln(𝑅 𝑒 )-1.41 1+1.32√ 𝜀 𝐷	and 𝐵 2 =	𝜀 3.7𝐷 𝑅 𝑒 + 2.51𝐵 1

Table 5 -

 5 2: Comparison of line flows in the electricity network -Case 5-A.

	Line	[95]	Electricity flow (kW + jkvar) Proposed
	i-j	(𝑃 𝑒𝑝 + 𝑗𝑃 𝑒𝑞 ) 𝑖𝑗	Losses	(𝑃 𝑒𝑝 + 𝑗𝑃 𝑒𝑞 ) 𝑖𝑗	Losses
	1 -2 26980.6+j15810.6 0.4352+j4351.8 26981+j15811	0.4352+j4351.8
	1 -3 27740.5+j11552.4 0310.2+j3101.9 23740+j11552	0310.2+j3101.9
	2 -3 -3454.6-j3541.2	12.7+j0127.4	-3454.6-j3541.2 12.7+j0127.4

Table 5 -

 5 3 : Comparison of hub parameters for heating network -Case 5-A.

	Hub	Nodal hydraulic head Heat generation and consumption (kW) Proposed [95] Proposed	Temperature ( °C) [95] Proposed
	i	𝐻 𝑠 (m)	𝐻 𝑟 (m)	𝐿 ℎ	𝑃 ℎ𝑔	𝐿 ℎ	𝑃 ℎ𝑔	𝑇 𝑠	𝑇 𝑟	𝑇 𝑠	𝑇 𝑟
	1	5517	10	0	28676.8	0	28647	120 48.6805	120 48.536
	2	10.666	5516.3	35000	0	35000	0	119.04	50 117.09	50
	3	4383.8	1143.2	20000	29000	20000	29000 123.541	49.534 123.54	49.04

Table 5 -

 5 4: Comparison of branch parameters in the heating network -Case 5-A. Proposed 𝒊 -𝒋 𝑚̇𝑖 𝑗 (kg/s) 𝑃ℎ 𝑖𝑗 (kW) Loss (kW) 𝑚̇𝑖 𝑗 (kg/s) 𝑃ℎ 𝑖𝑗 (kW Loss (kW)

	Branch		[95]	Heating network
	1 -2	64.6943 19175.967	890.1	65.962	19547	890.28
	1 -3	31.4533	9500.8	877	29.893	9100.2	873.09
	2 -3	-56.5335	-1674.1	909.7	-58.778	-16343	883.88

Table 5 -

 5 6: Input and outputs of the CHP plants at different hubs -Cases 5-B to 5-D.

	750
	650
	550
	450
	350
	250
	150
	50
	50

Table 5 -

 5 7: Hub-level operating parameters of the electricity and heating networks.

	(a) Case 5-B
	Electricity
	Hub

Table 5 -

 5 8: Selected branch-level operating parameters in the electricity and heating networks. The current per unit values are based on 4.16 kV line-to-line voltage and 500 kVA apparent power base values.

	(a) Case 5-B
	Hubs	Electricity

† 𝑇 𝑠(𝑖) ( o C) 𝑇 𝑠(𝑗) ( o C) 𝑇 𝑟(𝑖) ( o C) 𝑇 𝑟(𝑗) ( o C) *

Table 5 -

 5 10: Summary of the demand profile considered in Case 5-E.

			Electricity demand			Heat demand
	Hubs	Active lowest (kW) Peak (kW) Lowest (kvar) Peak (kvar) Reactive	lowest (kW) Peak (kW)	Load type
	2	8.68	67.10	5.38	41.6	2.61	77.09	Residential
	4	17.36	134.21	5.73	44.29	5.22	154.19	Residential
	6	50.68	199.50	24.33	95.76	2.00	54.64	Small hotel

Table 5 -

 5 11: Summary of the technologies considered in the two scenarios -Case 5-E.

			Technologies connected
	Hub 1	Hub 2	Hub 3	Hub 4	Hub 5	Hub 6
	Scenario I Slack hub HP2	CHP3	Boiler4, HP4 Solar PV HP6
	Scenario II Slack hub	---	CHP3, HP3(=HP2+HP4+HP6)	Boiler4	Solar PV	---

  -13, are quite the same. As there is no local reactive electricity generation, all the reactive power demands of the connected electric loads, the HPs and the CHP are covered by importing from the neighbourhood. The total reactive power demand of connected electric loads in both scenarios is 16.507 Mvarh. The reactive power consumption by the HPs is 2.431 Mvarh while the consumption by the CHP is 5.472 Mvarh. The reactive power lost in the network is about 290.394 kvarh in Scenario I and 272.877 kvarh in Scenario II. Thus, Scenario II is marginally more efficient than Scenario I from the electricity network point of view.

Table 6 -

 6 2: Optimal location and dispatches of the HPs and the HOB for Case 6-A.

	Scenarios	Energy carrier		HPs	HOB	CHP	PV	Imported Loss
		Location (hub numbers)	2	6	4	3	5	1	total
	I	Active electricity (kW) Reactive electricity (kvar) -9.42 19.45	24.67 -11.95 0.00 0.00	62.00 -30.03 0.00 50.00 350.38 239.90	17.45 25.71
		Heat (kW)	77.80	98.68	0.00	114.00 0.00	0.00	4.58
		Location (hub numbers)	2	2	6	3	5	1	total
	II	Active electricity (kW) Reactive electricity (kvar) -12.11 -9.36 25.00 19.33	0.00 0.00	62.00 -30.03 0.00 50.00 347.52 235.47	14.37 21.17
		Heat (kW)	100.00 77.32	0.00	114.00 0.00	0.00	5.38

Table 7 -

 7 1: Parameters describing a stratified thermal storage model

	Symbols	Definition
	𝐾 𝑠 𝑘 𝑤 𝑍 1 , 𝑍 2 and 𝑍 3 R storage 𝐴 𝑐 𝐴 𝑠1 , 𝐴 𝑠2 and 𝐴 𝑠3 𝑉 𝑠𝑡𝑜𝑟𝑎𝑔𝑒 𝑡 Δ𝑡 𝑇 1(𝑡-∆𝑡) , 𝑇 2(𝑡-∆𝑡) and 𝑇 3(𝑡-∆𝑡) 𝑇 1(𝑡) , 𝑇 2(𝑡) and 𝑇 3(𝑡) 𝑇 𝑎(𝑡) 𝑚 1 , 𝑚 2 and 𝑚 3 𝑚̇	Heat loss factor of the tank (W.m -2. K -1 ) Thermal conductivity of water (W.m -1 .K -1 ) Heights of layer 1, layer 2 and layer 3 (m) Radius (m) Cross-section area of the storage (m 2 ) Lateral surface areas of layer 1, layer 2 and layer 3 (m 2 ) Volume (m 3 ) Current time (s) Time step (s) Temperatures (K) of water in layer 1, layer 2 and layer 3 at time 𝑡 -Δ𝑡 Temperatures (K) of water in layer 1, layer 2 and layer 3 at time 𝑡 Ambient (soil) temperature (K) at time 𝑡 Mass (kg) of water in layer 1, layer 2 and layer 3 Mass flow rate (kg/s) going into or leaving from the tank

  |𝑉 𝑖 ||𝑉 𝑗 |(𝐺 𝑖𝑗 𝑐𝑜𝑠𝜃 𝑖𝑗 + 𝐵 𝑖𝑗 𝑠𝑖𝑛𝜃 𝑖𝑗 ), 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑗 ≠ 𝑖 (B-5)

	𝜕𝑃 𝑒𝑝 𝑖 𝜕𝜃 𝑗	= |𝑉 𝑖 ||𝑉 𝑗 |(𝐺 𝑖𝑗 𝑠𝑖𝑛𝜃 𝑖𝑗 -𝐵 𝑖𝑗 𝑐𝑜𝑠𝜃 𝑖𝑗 ),	𝑓𝑜𝑟 𝑎𝑙𝑙 𝑗 ≠ 𝑖		(B-1)
	𝜕𝑃 𝑒𝑝 𝑖 𝜕𝜃 𝑖	𝑁 = ∑ |𝑉 𝑖 | 𝑗=1,𝑗≠𝑖	|𝑉 𝑗 |(-𝐺 𝑖𝑗 𝑠𝑖𝑛𝜃 𝑖𝑗 + 𝐵 𝑖𝑗 𝑐𝑜𝑠𝜃 𝑖𝑗 ) = -𝑃 𝑒𝑞 𝑖	𝑖	-𝐵 𝑖𝑖 |𝑉 𝑖 | 2	(B-2)
	𝜕𝑃 𝑒𝑝 𝑖 𝜕|𝑉 𝑗 |	= |𝑉 𝑖 |(𝐺 𝑖𝑗 𝑐𝑜𝑠𝜃 𝑖𝑗 + 𝐵 𝑖𝑗 𝑠𝑖𝑛𝜃 𝑖𝑗 ),	𝑓𝑜𝑟 𝑎𝑙𝑙 𝑗 ≠ 𝑖			(B-3)
	𝜕𝑃 𝑒𝑝 𝑖 𝜕|𝑉 𝑖 |	= ∑ |𝑉 𝑗 | 𝑁 𝑗=1,𝑗≠𝑖	(𝐺 𝑖𝑗 𝑐𝑜𝑠𝜃 𝑖𝑗 + 𝐵 𝑖𝑗 𝑠𝑖𝑛𝜃 𝑖𝑗 ) + 2|𝑉 𝑖 |𝐺 𝑖𝑖 =	𝑃 𝑒𝑝 𝑖 |𝑉 𝑖 |	+ 𝐺 𝑖𝑖 |𝑉 𝑖 |	(B-4)
	𝜕𝑃 𝑒𝑞 𝑖 𝜕𝜃 𝑗 = -𝜕𝑃 𝑒𝑞 𝑖 𝑗=1,𝑗≠𝑖 𝑁 𝜕𝜃 𝑖 = ∑ |𝑉 𝑖 |	|𝑉 𝑗 |(𝐺 𝑖𝑗 𝑐𝑜𝑠𝜃 𝑖𝑗 + 𝐵 𝑖𝑗 𝑠𝑖𝑛𝜃 𝑖𝑗 ) = 𝑃 𝑒𝑝 𝑖 -𝐺 𝑖𝑖 |𝑉 𝑖 | 2	(B-6)
	𝜕𝑃 𝑒𝑞 𝑖 𝜕|𝑉 𝑗 |	= |𝑉 𝑖 |(𝐺 𝑖𝑗 𝑠𝑖𝑛𝜃 𝑖𝑗 -𝐵 𝑖𝑗 𝑐𝑜𝑠𝜃 𝑖𝑗 ) ,	𝑓𝑜𝑟 𝑎𝑙𝑙 𝑗 ≠ 𝑖			(B-7)
	𝜕𝑃 𝑒𝑞 𝑖 𝜕|𝑉 𝑖 |	= ∑ |𝑉 𝑗 | 𝑁 𝑗=1,𝑗≠𝑖	(𝐺 𝑖𝑗 𝑠𝑖𝑛𝜃 𝑖𝑗 -𝐵 𝑖𝑗 𝑐𝑜𝑠𝜃 𝑖𝑗 ) -2|𝑉 𝑖 |𝐵 𝑖𝑖 =	𝑃 𝑒𝑞 𝑖 |𝑉 𝑖 |	-𝐵 𝑖𝑖 |𝑉 𝑖 | 2	(B-8)

: Equations for the first row of the Jacobian matrix:

				𝜕(∆𝑃 𝑒𝑝 ) 𝑖 𝜕𝑇 𝑟(𝑗)	= 𝐶 𝑒𝑝(ℎ) 𝑖	𝜕𝑃 ℎ 𝑖 𝜕𝑇 𝑟(𝑗)	(B-20)
								B.2:
				𝝏(∆𝑷 𝒆𝒑 )			𝝏(∆𝑷 𝒆𝒑 )	𝝏(∆𝑷 𝒆𝒑 )	𝝏(∆𝑷 𝒆𝒑 )	𝝏(∆𝑷 𝒆𝒑 )	𝝏(∆𝑷 𝒆𝒑 )
				𝝏𝜽 𝝏(∆𝑷 𝒆𝒒 )			𝝏|𝑽| 𝝏(∆𝑷 𝒆𝒒 )	𝝏𝒎 ̇𝒑𝒊𝒑𝒆 𝝏(∆𝑷 𝒆𝒒 )	𝝏𝑯 𝝏(∆𝑷 𝒆𝒒 )	𝝏𝑻 𝒔 𝝏(∆𝑷 𝒆𝒒 )	𝝏𝑻 𝒓 𝝏(∆𝑷 𝒆𝒒 )
	[ ∆𝑯 𝒑𝒊𝒑𝒆 ∆𝑷 𝒆𝒑 ∆𝑷 𝒆𝒒 ∆𝑷 𝒉 ∆𝑬 𝒔 ∆𝑬 𝒓 ]	=	𝝏𝜽 𝝏(∆𝑷 𝒉 ) 𝝏𝜽 𝝏(∆𝑯 𝒑𝒊𝒑𝒆 ) 𝝏𝜽 𝝏(∆𝑬 𝒔 )			𝝏|𝑽| 𝝏(∆𝑷 𝒉 ) 𝝏|𝑽| 𝝏(∆𝑯 𝒑𝒊𝒑𝒆 ) 𝝏|𝑽| 𝝏(∆𝑬 𝒔 )	𝝏𝒎 ̇𝒑𝒊𝒑𝒆 𝝏(∆𝑷 𝒉 ) 𝝏𝒎 ̇𝒑𝒊𝒑𝒆 𝝏(∆𝑯 𝒑𝒊𝒑𝒆 ) 𝝏𝒎 ̇𝒑𝒊𝒑𝒆 𝝏(∆𝑬 𝒔 )	𝝏𝑯 𝝏(∆𝑷 𝒉 ) 𝝏𝑯 𝝏(∆𝑯 𝒑𝒊𝒑𝒆 ) 𝝏𝑯 𝝏(∆𝑬 𝒔 )	𝝏𝑻 𝒔 𝝏(∆𝑷 𝒉 ) 𝝏𝑻 𝒔 𝝏(∆𝑯 𝒑𝒊𝒑𝒆 ) 𝝏𝑻 𝒔 𝝏(∆𝑬 𝒔 )	𝝏𝑻 𝒓 𝝏(∆𝑷 𝒉 ) 𝝏𝑻 𝒓 𝝏(∆𝑯 𝒑𝒊𝒑𝒆 ) 𝝏𝑻 𝒓 𝝏(∆𝑬 𝒔 )	[	∆𝜽 ∆|𝑽| ∆𝒎 ̇𝒑𝒊𝒑𝒆 ∆𝑻 𝒓 ] ∆𝑻 𝒔 ∆𝑯	(B-14)
				𝝏𝜽 𝝏(∆𝑬 𝒓 )			𝝏|𝑽| 𝝏(∆𝑬 𝒓 )	𝝏𝒎 ̇𝒑𝒊𝒑𝒆 𝝏(∆𝑬 𝒓 )	𝝏𝑯 𝝏(∆𝑬 𝒓 )	𝝏𝑻 𝒔 𝝏(∆𝑬 𝒓 )	𝝏𝑻 𝒓 𝝏(∆𝑬 𝒓 )
			[	𝝏𝜽			𝝏|𝑽|	𝝏𝒎 ̇𝒑𝒊𝒑𝒆	𝝏𝑯	𝝏𝑻 𝒔	𝝏𝑻 𝒓 ]
	𝜕(∆𝑃 𝑒𝑝 ) 𝑖 𝜕𝜃 𝑗	= 𝐶 𝑒𝑝(𝑒𝑝) 𝑖	𝜕𝑃 𝑒𝑝 𝑖 𝜕𝜃 𝑗	+ 𝐶 𝑒𝑝(𝑒𝑞) 𝑖	𝜕𝑃 𝑒𝑞 𝑖 𝜕𝜃 𝑗	(B-15)
	𝜕(∆𝑃 𝑒𝑝 ) 𝑖 𝜕|𝑉| 𝑗	= 𝐶 𝑒𝑝(𝑒𝑝) 𝑖	𝜕𝑃 𝑒𝑝 𝑖 𝜕|𝑉| 𝑗	+ 𝐶 𝑒𝑝(𝑒𝑞) 𝑖	𝜕𝑃 𝑒𝑞 𝑖 𝜕|𝑉| 𝑗	(B-16)
	𝜕(∆𝑃 𝑒𝑝 ) 𝑖 𝜕𝑚̇𝑝 𝑖𝑝𝑒(𝑚)	= 𝐶 𝑒𝑝(ℎ) 𝑖	𝜕𝑃 ℎ 𝑖 𝜕𝑚̇𝑝 𝑖𝑝𝑒(𝑚)	(B-17)
	𝜕(∆𝑃 𝑒𝑝 ) 𝑖 𝜕𝐻 𝑗	= 0								(B-18)
	𝜕(∆𝑃 𝑒𝑝 ) 𝑖 𝜕𝑇 𝑠(𝑗)	= 𝐶 𝑒𝑝(ℎ) 𝑖	𝜕𝑃 ℎ 𝑖 𝜕𝑇 𝑠(𝑗)				(B-19)

B.

[START_REF] Joyeux | Household energy consumption versus income and relative standard of living: A panel approach[END_REF]

Equations for the second row of the Jacobian matrix:

	𝜕(∆𝑃 𝑒𝑞 ) 𝑖 𝜕𝜃 𝑗	= 𝐶 𝑒𝑞(𝑒𝑝) 𝑖	𝜕𝑃 𝑒𝑝 𝑖 𝜕𝜃 𝑗	+ 𝐶 𝑒𝑞(𝑒𝑞) 𝑖	𝜕𝑃 𝑒𝑞 𝑖 𝜕𝜃 𝑗	(B-21)
	𝜕(∆𝑃 𝑒𝑞 ) 𝑖 𝜕|𝑉| 𝑗	= 𝐶 𝑒𝑞(𝑒𝑝) 𝑖	𝜕𝑃 𝑒𝑝 𝑖 𝜕|𝑉| 𝑗	+ 𝐶 𝑒𝑞(𝑒𝑞) 𝑖	𝜕𝑃 𝑒𝑞 𝑖 𝜕|𝑉| 𝑗	(B-22)
	𝜕(∆𝑃 𝑒𝑞 ) 𝑖 𝜕𝑚̇𝑝 𝑖𝑝𝑒(𝑚)	= 𝐶 𝑒𝑞(ℎ) 𝑖	𝜕𝑃 ℎ 𝑖 𝜕𝑚̇𝑝 𝑖𝑝𝑒(𝑚)	(B-23)
	𝜕(∆𝑃 𝑒𝑞 ) 𝑖 𝜕𝐻 𝑗	= 0				(B-23)
	𝜕(∆𝑃 𝑒𝑞 ) 𝑖 𝜕𝑇 𝑠(𝑗)	= 𝐶 𝑒𝑞(ℎ) 𝑖	𝜕𝑃 ℎ 𝑖 𝜕𝑇 𝑠(𝑗)	(B-25)
	𝜕(∆𝑃 𝑒𝑞 ) 𝑖 𝜕𝑇 𝑟(𝑗)	= 𝐶 𝑒𝑞(ℎ) 𝑖	𝜕𝑃 ℎ 𝑖 𝜕𝑇 𝑟(𝑗)	(B-26)

B.

3: Equations for the third row of the Jacobian matrix:

	𝜕(∆𝑃 ℎ ) 𝑖 𝜕𝜃 𝑗	= 𝐶 ℎ(𝑒𝑝) 𝑖	𝜕𝑃 𝑒𝑝 𝑖 𝜕𝜃 𝑗	+ 𝐶 ℎ(𝑒𝑞) 𝑖	𝜕𝑃 𝑒𝑞 𝑖 𝜕𝜃 𝑗	(B-27)
	𝜕(∆𝑃 ℎ ) 𝑖 𝜕|𝑉| 𝑗	= 𝐶 ℎ(𝑒𝑝) 𝑖	𝜕𝑃 𝑒𝑝 𝑖 𝜕|𝑉| 𝑗	+ 𝐶 ℎ(𝑒𝑞) 𝑖	𝜕𝑃 𝑒𝑞 𝑖 𝜕|𝑉| 𝑗	(B-28)
	𝜕(∆𝑃 ℎ ) 𝑖 𝜕𝑚̇𝑝 𝑖𝑝𝑒(𝑚)	= 𝐶 ℎ(ℎ) 𝑖	𝜕𝑃 ℎ 𝑖 𝜕𝑚̇𝑝 𝑖𝑝𝑒(𝑚)	(B-29)
	𝜕(∆𝑃 ℎ ) 𝑖 𝜕𝐻 𝑗	= 0				(B-30)
	𝜕(∆𝑃 ℎ ) 𝑖 𝜕𝑇 𝑠(𝑗)	= 𝐶 ℎ(ℎ) 𝑖	𝜕𝑃 ℎ 𝑖 𝜕𝑇 𝑠(𝑗)	(B-31)
	𝜕(∆𝑃 ℎ ) 𝑖 𝜕𝑇 𝑟(𝑗)	= 𝐶 ℎ(ℎ) 𝑖	𝜕𝑃 ℎ 𝑖 𝜕𝑇 𝑟(𝑗)	(B-32)

Table C -

 C 2: Carrier pipe, insulation and soil parameters for the heating networks -Cases 5-B to 5-D

	Hubs				Carrier Pipe Parameters		Insulation	Outer Jacket		Soil
	from to	L (m)	𝑅 1 (mm)	𝑅 2 (mm)	𝑘 2 (W/mK)	e (mm)	𝑅 3 (mm)	k3 (W/mK)	𝑅 4 (mm)	𝑘 4 (W/mK)	Z (m)	𝑘 𝑠 (W/mK)	flow (kg/s)
	1	2	3200	26.95	30.15	40.5	0.05	67	0.027	70	0.40	0.5	1.15	10 *
	2	4	2500	26.95	30.15	35.6	0.06	67	0.027	70	0.40	0.5	1.15	10 *
	3	2	1300	26.95	30.15	45.1	0.09	67	0.027	70	0.40	0.5	1.15	10 *
	3	4	1800	26.95	30.15	42.0	0.10	67	0.027	70	0.40	0.5	1.15	10 *
	1	5	3800	26.95	30.15	38.0	0.08	67	0.027	70	0.40	0.5	1.15	10 *
	4	5	2000	26.95	30.15	40.0	0.10	67	0.027	70	0.40	0.5	1.15	10

* Table C-3: Generation and demands of electricity, heat and fuel power at different hubs -Cases 5-B to 5-D Hub Electricity

network hub data District heating network hub data Fuel Type

  V (pu) 𝜃 ( 0 ) 𝑃 𝑒𝑝𝑔 (kW) 𝑃 𝑒𝑞𝑔 (kvar) 𝐿 𝑒𝑝 kW) 𝐿 𝑒𝑞 (kvar) Type 𝑇 𝑠 ( 0 C) 𝑇 𝑟 ( 0 C) 𝑃 ℎ𝑔 (kW) 𝐿 ℎ (kW) 𝐻 𝑠 (m) 𝐻 𝑟 (m) 𝑃 𝑓𝑔 (kW) For electricity Type 1 is slack, 2 is PV and 3 is PQ while for DHN type 1 is slack, 2 is temp_supply and 3 is temp_return Equation (C-1) represents the general format of the coupling matrix and specific coupling matrixes for each hub are given below.𝐶 𝑒𝑝(𝑒𝑝) 𝐶 𝑒𝑝(𝑒𝑞) 𝐶 𝑒𝑝(ℎ) 𝐶 𝑒𝑝(𝑓) 𝐶 𝑒𝑞(𝑒𝑝) 𝐶 𝑒𝑞(𝑒𝑞) 𝐶 𝑒𝑞(ℎ) 𝐶 𝑒𝑞(𝑓) 𝐶 ℎ(𝑒𝑝) 𝐶 ℎ(𝑒𝑞) 𝐶 ℎ(ℎ) 𝐶 ℎ(𝑓)

	1 1	1	0 0 *	0 *	25	0	1	80	20 *	0 *	30	135 1	0 *
	2 3	1 *	0 * 0	0	160	60	3 † 65 *	30	50	250	50 *	1 *	50
	3 2	1.05 0 * 50	20 *	200	130	2 † 90	25 *	35	50	50 *	1 *	500
	4 3	1 *	0 * 0	0	50	80	3	70 *	35	20	220	50 *	1 *	50
	5 3	1 *	0 * 0	0	90	40	3	100 * 35	0	150	50 *	1 *	55
	NB: 𝐶 𝑔𝑒𝑛𝑒𝑟𝑎𝑙 =

* A line-to-line voltage base value of 4.16 kV is considered in the simulation.

* Most of currently operating DHNs (categorized as 3 rd generations) supply heat using pressurised water in the range of 70 -100 °C while the future 4 th generation are expected to use supply temperatures as low as 50 -60 °C[START_REF] Lund | 4th Generation District Heating (4GDH): Integrating smart thermal grids into future sustainable energy systems[END_REF]. The corresponding return temperatures are lower than 45 °C and 25 °C, respectively. In this thesis, temperatures in the ranges of both 3 rd and 4 th generation DHNs are considered.

Acknowledgements | xxv

Acknowledgements

List of publications

Six peer-reviewed scientific publications are made during the course of this PhD program. All of them are original research papers. Four of them are either published on or submitted to international journals while the other two are presented at international conferences. Set maximum iterations for outer PSO as MaxIter1 and start with Iter1=1, and particle i=1 Set maximum iteration for inner PSO as MaxIter2 and start with Iter2=1, and particle j=1 Fig. 6-1. Flow chart of the nested PSO applied for optimal placement and sizing of HPs and a HOB -Cases 6-A and 6-B.

Hubs 3 and 4 are sources of heat as they are indicated by positive nodal mass flows in Table 6-12 while hubs 2 and 5 are consumers. The return temperature of both hubs 2 and 5 are at 30 °C which is the lower limit. It implies that the lower the return temperature from the consumers the more economical the DHN will be.

The supply temperatures of hubs 3 and 4 are at 95 °C and 67.76 °C, respectively. Should the cost of pumping energy be not significant, both supply temperatures were expected to be equal to the minimum limit of 60 °C. However, both of them are higher than the lower limit and hub 3 is rather at the maximum limit. This shows that decreasing the sources temperatures further leads to not only a decrease in distribution loss, but also an increase in a mass flow rate, which in turn requires higher pumping energy. The optimisation found that the cost of the additional pumping energy is more expensive than the cost of heat loss that could be saved by decreasing the supply temperatures further than those values for this particular case study.

As it can be calculated from Newton-Raphson numerical method is implemented in this thesis aiming to show the more comprehensive approach of solving the integrated load flow problem. Other numerical techniques, such as Gauss-Seidel and Quasi-Newton methods, can be tested so that the computational cost can be reduced while simulating larger network sizes. Algorithms that are specific to radial distribution networks can also be investigated using the Newton-Raphson method as a benchmark. Similarly, the simple greedy algorithm used for topology optimisation can be further tuned to reduce the likelihood of sub-optimal solutions. It can also be extended to larger case studies of meshed electricity and heat distribution networks.

Energy demands at the end users are assumed to be very accurate forecasts. In addition, intermittent renewable energy sources, such as solar PV and wind power generations are included in different case studies. The HOMER [START_REF] Homer Energy | HOMER -Hybrid Renewable and Distributed Generation System Design Software[END_REF] simulation tool is used to estimate their generation outputs. The forecasting procedure, however, is not discussed in this thesis. These assumptions can be improved by considering uncertainties, which is another research area where the proposed tool can be used as a base. The tool can also be used to conduct sensitivity analysis on MES as a function of renewables' availability, variations in price signals and location of coupling technologies (e.g. HP vs CHP as a function of electricity and fuel price variations).

The other limitation is on the assumption that the control system in both heating and electricity distribution networks are capable of realizing the load flow and thermo-economic optimisation results. For example, the circulation pumps and pressure valves of the heating network are not explicitly included in the model. Also the control system on the secondary side of each heat substation is assumed to have the capacity to handle the temperature and mass flow rates considered in the primary side. Any deviations from the reality can be addressed by adjusting the constraints of the operating parameters of the distribution network, such as temperature and mass flow rates.

The distribution networks and the associated energy technologies considered in all case studies are assumed operational. As a result, the operating cost is considered as a main variable in the objective functions of the optimisation case studies. This can be complemented by considering the emission and investment cost of the network and coupling technologies which can be then used for planning future energy systems or expansion of existing ones.

Yet another limitation is with regard to the thermal inertia in the heat distribution networks. The model presented in this thesis considers hourly steady state operation of both heating and electricity networks. Though the electricity network reaches to steady state in fraction of seconds, the thermal networks may take from minutes up to hours to reach a steady state depending on their sizes. In this thesis, the case studies considered are small enough to assume hourly steady state operation and neglect the heating network inertia. For very large heating networks, however, consideration of thermal inertia gives more accurate results from operational control point of view. In that regard, the way the stead state models of the corresponding networks are integrated can be used as a base to integrate the dynamic models of the electricity and heat distribution networks. That opens a wider research area in the integration of smart grids and smart thermal networks.

Appendices Appendix A. Derivation of admittance matrix and electrical power injection equations

In per unit system (the same equation for single phase and three phase system), the apparent power injected at bus 𝑖 (𝑆 𝑒(𝑖) ) in complex number notation is defined by equation (A-1) [START_REF] Glover | Power System Analysis & Design, SI Version[END_REF].

where 𝑉 𝑖 is the voltage at bus 𝑖 while 𝐼 𝑖 * is the conjugate of the current injected at bus 𝑖. According to Kirchhoff's current law, the current injected at bus 𝑖 is the sum of all outgoing current flows from the same bus. Each outgoing current flow is then expressed in terms of Ohm's law as a function of line admittance (𝑦 𝑖𝑗 ) and voltage difference at the two ends of the line (𝑉 𝑖 -𝑉 𝑗 ) as shown in equation (A-2).

The term 𝑦 𝑖𝑖 in equation (A-2) refers to half of the susceptance of the line aggregated at the two ends of the line (usually designated by 𝑏/2), and 𝑁 is the total number of buses in the network. Equation (A-2) can be further rearranged as shown in equations (A-3). and 𝑌 𝑖𝑗 = -𝑦 𝑖𝑗 In a more compact form, the current injection described in equation (A-4) can be rewritten as shown in equation (A-5) where 𝑰, 𝒀 and 𝑽 are matrices representing current injections at all buses, network admittance matrix and voltages at all buses.

The elements of network admittance matrix can be represented in rectangular form (𝑌 𝑖𝑗 = 𝐺 𝑖𝑗 + 𝑗𝐵 𝑖𝑗 ) while the voltages are represented in polar form (𝑉 𝑖 = |𝑉 𝑖 |∠𝜃 𝑖 ). Hence, the per-unit apparent/complex power can be formulated by substituting equation (A-5) into (A-1). 𝑃 𝑒𝑝(𝑖) and 𝑃 𝑒𝑞(𝑖) are the active and reactive components of injected power at bus 𝑖. After substitution for 𝑉 and 𝑌 and separating real and imaginary parts, the active power and reactive power injections at bus 𝑖 are given by equations (A-7) and (A-8), respectively.

𝜃 𝑖𝑗 is the voltage angle difference between nodes 𝑖 and 𝑗 (𝜃 𝑖 -𝜃 𝑗 ).

For a transmission line with a tap changing transformer, the admittance matrix needs to be updated by considering the pi-equivalent of the transformer.

B.4: Equations for the fourth row of the Jacobian matrix:

B.5: Equations for the fifth row of the Jacobian matrix:

In equations (B-41), if the flow in pipe 𝑛 is positive and getting into node 𝑖 equation (B-41a) is used. Instead, if the pipe flow is negative (opposite to the positive direction) and is still flowing into node 𝑖, then equation (B-41b) is used. Otherwise, equation (B-41c) is applied.