Laurent Lefevre

Robin Roche

Mr Than

Hung Pham

This thesis is the result of a three-year work within the laboratory LCIS of Grenoble INP based in Valence, France.

Résumé

L'objectif principal de la thèse est la résolution du problème de gestion d'énergie dans un microréseau maillé en courant continu (CC). Le microréseau à l'étude est composé d'une source d'énergie renouvelable (panneaux solaires), d'un système de stockage d'énergie, lui-même composé de batteries au plomb-acide, de convertisseurs CC et de charges (appareils domestiques ou professionnels, véhicules électriques, etc.). Le microréseau est connecté au réseau externe qui, avec les panneaux solaires, constituent les sources principales du système. Le microréseau en CC peut être déconnecté du réseau externe et fonctionner ainsi indépendamment en évitant les coupures de courant occasionnelles. Le réseau de transmission central, appelé bus CC commun, est un réseau de haute tension. Les pertes de puissance dans les lignes de transmission sont prises en compte. L'équilibrage de la puissance, la minimisation des coûts, la dissipation de la puissance, la planification de la batterie, le contrôle des interrupteurs dans les convertisseurs seront pris en compte.

Les microréseaux sont des systèmes physiques compliqués et de nombreux problèmes se posent et restent ouverts. Ces problèmes concernent notamment la distribution d'énergie, la qualité d'énergie, la topologie, le dimensionnement des sources d'énergie renouvelables et des éléments de stockage, le type des composants et la structure. La dynamique est complexe avec des échelles temporelles multiples, des non linéarités, des contraintes, des problèmes de stabilité et des incertitudes. Par conséquent, la mise en oeuvre d'une stratégie de commande optimale est difficile et des approches de modélisation et de commande spécifiques doivent être considérées.

L'approche hamiltonienne à port est une méthode de modélisation qui décrit la conservation de la puissance et de l'énergie entre les éléments physiques du réseau. La représentation d'état sous forme hamiltonienne peut être générée à partir du Bond graph associé, qui est une représentation graphique intuitive de type circuit du système consideré. Un modèle hamiltonien à port est proposé pour le micro-réseau en CC considéré, comprenant les représentations hamiltoniennes à port pour chaque composant du système séparément. Le modèle hamiltonien à port est inversé en utilisant une méthode de platitude différentielle qui conduit à la représentation plate de chaque partie du micro-réseau (PV, ES, charges, tous connectés aux convertisseurs CC). Avec la dynamique inverse, les états et les entrées de commande sont dérivés en fonction des sorties plates choisies et sont utilisés par la suite dans la formulation du problème d'optimisation. La dynamique multi-échelle est traitée par une approche de commande hiérarchique divisée en trois niveaux. Tout d'abord, au premier niveau, des profils optimaux sont générés pour une approche d'optimisation qui vise à l'équilibrage de la puissance, la réduction de coût et la diminution de la dissipation. La paramétrisation B-spline est considérée pour la validation des contraintes à temps continu dans ce problème d'optimization et pour la génération de trajectoires plates. Au niveau secondaire, un contrôleur de type Model Predictive Control est développé pour suivre les profils optimaux obtenus dans le niveau supérieur. Des perturbations sont prises en compte et les erreurs du système sont analysées. Enfin, au niveau bas, les profils de suivi sont considérés comme consignes pour la régulation des convertisseurs CC. L'approche proposée est validée à travers des simulations approfondies sur un microréseau maillé en CC utilisant des profils réels pour les charges, la température extérieure et le rayonnement solaire.

Notations

A description of some notions is provided in order to facilitate the coherence and the understanding of the techniques discussed in the text.

Bond graph [START_REF] Karnopp | System dynamics: modeling, simulation, and control of mechatronic systems[END_REF]: It is a graph-oriented tool which describes the dynamics of complex physical systems. It allows the conversion from a graphical to a state-space representation maintaining the power conservation among the elements.

Depth of discharge of a battery: It indicates the percentage of the capacity removed with respect to the overall capacity of the battery [START_REF] Spiers | Batteries in pv systems[END_REF].

State of charge of a battery: It indicates the percentage of the capacity which remains in the battery with respect to the overall capacity [START_REF] Spiers | Batteries in pv systems[END_REF].

Duty cycle: It describes the time of an active signal over one period.

Fill factor: It is defined by the difference between the maximum power and the theoretical expected power generated by a solar cell. Down-conversion in a DC/DC converter: It is the process during which the converter produces an output voltage lower than the input voltage.

Up-conversion in a DC/DC converter: It is the process during which the converter produces an output voltage higher than the input voltage.

Jumper: It is a small electronic device which operates as a conductor and allows the electricity to flow trough specific paths of the electrical circuits. It functions, also, as a switch which can close or open the path.

Manifold: It is a topological space which is locally similar to an Euclidean space. A n-manifold M is a topological space locally homeomorphic with the Euclidean space at the n point [START_REF] Cannon | The recognition problem: what is a topological manifold?[END_REF].

Tangent space: Consider a n-manifold M and x a point of M . The R n tangents attach to the x point form the tangent space of M , which is written as T x M . If γ is a smooth curve which passes from x, then its derivative is a vector in T x M (Weisstein, 2019).

Tangent bundle: It is the manifold T M , which gathers all the tangent vectors in T x M for every point x of the n-manifold M (Berger and Gostiaux, 2012).

One-form: An one-form, ω, is a mapping of the tangent bundle T M of the manifold M , ω : T M -→ R. Locally, it is a linear map of the tangent space T x M , ω x : T x M -→ R, described by a linear combination of differentials of the points x: ω x = Nomenclature Bond graph (causal and bicausal) and port-Hamiltonian systems

Introduction

The demand for independent power systems connected to the main grid has been continuously increasing over the past few years. Their existence plays a vital role since it leads to a more efficient operation for the main grid. They minimize the power losses, store the energy excess and manage load sharing. These systems include distributed generators (DG), renewable resources (e.g. solar panels (PV) or wind turbines) and storage systems (batteries, ultra-capacitors), linked to a variety of loads (industrial or residential) [START_REF] Planas | Ac and dc technology in microgrids: A review[END_REF]. In this regard, microgrids (Fig. 1.0.1) have gained important attention and they are considered as suitable systems able to contain DGs within the existing electrical grid system. Their appearance ensures a reliable power distribution and, hence, their existence into the grid operation support have become necessary [START_REF] Lotfi | Hybrid ac/dc microgrid planning[END_REF]. Although the concept of microgrids seems to have become increasingly popular in the latest years, their applications appeared long ago. In 1882, Tomas Edison opened the first commercial central power plant in Manhattan, USA, a self-maintained electrical system powered by coal with an energy storage (ES) unit and a limited distribution DC network [START_REF] Escobar | A hamiltonian viewpoint in the modeling of switching power converters[END_REF][START_REF] Hughes | Networks of power: electrification in Western society, 1880-1930[END_REF]. Afterwards, in the 90's, several projects in USA and in Europe held by the CERTS (Consortium for Electric Reliability Technology Solutions) [START_REF] Lasseter | Integration of distributed energy resources[END_REF] and the CIGRE (International Council of Large Electric Systems) [START_REF] Hatziargyriou | Microgrids-large scale integration of microgeneration to low voltage grids[END_REF] implemented some of the first modern grid-connected microgrids. This concept can incorporate multiple distributed energy resources and, in this way, create autonomous systems working in islanded mode under different control strategies [START_REF] Hirsch | Microgrids: A review of technologies, key drivers, and outstanding issues[END_REF]. Since then, several examples of microgrids can be found all 1 around the world, like the Kythnos Island in Greece, the MVV 1 Residential demonstration in Germany and the Continuon MV/LV2 microgrid in Netherlands [START_REF] Katiraei | Microgrids management[END_REF]Moehrke and Myrzik, 2014). Since the microgrids are complex power systems composed by a large number of sources, a huge amount of power electronic converters, a set of numerous loads and a multi-line transmission network, many issues arise for further study.

An important aspect is the type of the transmission network: AC or DC (Fig. 1.0.2). Despite the domination of AC transmission networks, the interest on DC microgrids grew lately as a result of the constant development and production of the DC equipment for renewable sources (e.g. PVs), storage systems (e.g. lead-acid or lithium-ion batteries) and loads (e.g. electrical vehicles, elevator and escalator systems and various smaller DC loads such as computers, LED lights and the like). Several works have shown that DC microgrids may have numerous advantages over an AC distribution. For instance in [START_REF] Whaite | Power quality in dc power distribution systems and microgrids[END_REF] they mention that a properly structured DC microgrid provides better power quality with ES devices and DC loads integration without frequency or phase control. Furthermore, easier integration of DC power sources and compact power electronics equipment are some of the other advantages according to [START_REF] Jin | Next-generation shipboard dc power system: Introduction smart grid and dc microgrid technologies into maritime electrical netowrks[END_REF]. Besides, the power losses minimization obtained by reducing the number of converters, avoidance of the multiple DC buses synchronization and an effective power delivery towards DC loads, such as electrical vehicles (EVs), are other beneficial aspects as reported by Lotfi and Khodaei (2017a) and [START_REF] Zubieta | Are microgrids the future of energy?: Dc microgrids from concept to demonstration to deployment[END_REF]. In general, the main goal is always to handle the load balancing problem [START_REF] Kanchev | Energy management and operational planning of a microgrid with a pv-based active generator for smart grid applications[END_REF] minimizing the cost, reducing the power losses and preventing the fault occurrences (e.g. power outages and lines under faults). However, the reliability, the stability and the quality of the power distribution in microgrids continue to be critical issues for the researcher community. Microgrids, as a group of distributed energy resources (DERs), interconnected loads and a large amount of transmission lines, are complex dynamical systems and up to date numerous problems have arisen and remain open for investigation: Firstly, an important factor is the topology of the system. There are two main types of topologies: the radial and the meshed topology. In the radial topology, also called point-to-point component interconnection, the electricity is distributed from an energy source to a number of loads. The disadvantage of this topology is that as the electricity passes through the power lines, its strength ends up reduced at the side of the load [START_REF] Flick | Securing the Smart Grid: Next Generation Power Grid Security[END_REF]. Additionally, in the case of power outages, the whole system gets influenced and suffers from constant power interruptions. Contrariwise, a meshed topology offers greater reliability since the electricity can be transmitted through a variety of sources and transmission lines. Thus, the efficiency and the continuity of the power transmission can be guaranteed and the safe operation of the system can be ensured [START_REF] Flick | Securing the Smart Grid: Next Generation Power Grid Security[END_REF][START_REF] Guerrero | Advanced control architectures for intelligent microgrids, part i: decentralized and hierarchical control[END_REF][START_REF] Haileselassie | Power system security in a meshed north sea hvdc grid[END_REF].

Another important aspect is the sizing of the renewable sources and the energy storing elements. Sizing algorithms have been developed which consider optimal power sizing of the PVs or the wind turbines depending on the economical cost, on the fuel emission minimization, the ES system or the type of the DC loads, like EVs [START_REF] Atia | Sizing and analysis of renewable energy and battery systems in residential microgrids[END_REF][START_REF] Iovine | Nonlinear control of a dc microgrid for the integration of photovoltaic panels[END_REF][START_REF] Kumar | Renewable energy based microgrid system sizing and energy management for green buildings[END_REF]. Besides, the ES system plays an important role in microgrids and ensures the efficient operation of the system [START_REF] Fossati | A method for optimal sizing energy storage systems for microgrids[END_REF][START_REF] Khorramdel | Optimal battery sizing in microgrids using probabilistic unit commitment[END_REF]. Its proper sizing can provide a backup energy source in case of power transfer interruption from the utility grid (UG) with concurrent inability of the renewable sources system to support the consumers' demand.

Another critical issue is the configuration of the energy management system (EMS), which manages the power balancing problem, minimizes the power dissipation and ameliorates the power quality. The EMS deals with the appropriate use of the energy sources and decides the proper power flow distribution to satisfy the consumers' demand [START_REF] Shayeghi | A survey on microgrid energy management considering flexible energy sources[END_REF]). An EMS can operate under a centralized or a decentralized strategy framework [START_REF] Meng | Microgrid supervisory controllers and energy management systems: A literature review[END_REF]. Under a centralized supervision, the system takes global optimal decisions collecting information from the components of the microgrid. Contrariwise, in the decentralized management, local supervision is possible which allows flexibility during operation and isolation of the events under fault. Furthermore, off-line and on-line EMS approaches have been developed which either predict the functional behavior of the system or study the real-time operation of the microgrid respectively [START_REF] Shi | Real-time energy management in microgrids[END_REF].

To deal with the above enumerated problems, a modeling and a supervision approach are necessary depending on the architecture of the system and the control objectives. In the following, an overview for the hierarchical control approach will be presented focusing mostly on the energy management of the meshed DC microgrids.

Hierarchical control

In microgrids, many matters emerge concerning multiple dynamics, nonlinearities, constraints, reliability, stability and so on. To cope with all these challenges, hierarchical control strategies have been developed which handle the aforementioned problems by decomposing them into smaller subproblems. the second level (secondary level), which is responsible for the voltage regulation alleviating the disturbances. Besides, in the case of multiple DC buses, this controller synchronizes also the connection among them; the first level (primary control), which deals with the local supervision of the power converters linked to the sources and the loads. This controller handles also the voltage and current sharing among them.

Since this approach can manage different control objectives at the same time, it appears to be a suitable and advantageous tool for DC microgrid control. Some of the benefits are enlisted below: handling of both fast and slow dynamics of the system. The slow dynamics accounts for the sources and the power flow. On the other hand, the fast dynamics describes the switching activity in the DC/DC converters; integration of different control methods at every level to deal with the nonlinear dynamics under constraints. For instance, at one level, the optimal power profiles can be generated and, at another level, the error dynamics can be studied for discrepancies' mitigation [START_REF] Pham | Power balancing in a dc microgrid elevator system through constrained optimization[END_REF]Zafeiratou et al., 2018a); considering energy management strategies to satisfy economical and technological criteria, such as minimization of the power losses, amelioration of the microgrid's architecture and the like; study of each component separately which gives the possibility of local supervision without influencing the rest of the system (Cairoli and Dougal, 2017); taking into account stability and robustness of the system, which guarantee the continuous and secure operation [START_REF] Guo | Stability analysis and damping enhancement based on frequency-dependent virtual impedance for dc microgrids[END_REF][START_REF] Herrera | Stability analysis and controller design of dc microgrids with constant power loads[END_REF]; reconfiguration of the system under unexpected faulted events, when, for example, a line stops transmitting power or a source is not capable of delivering the requested amount of energy.

Certainly, the advantages of the DC microgrids facilitate the implementation of the hierarchical control, since synchronization or reactive power issues no longer exist. Hereinafter, some significant factors will be described which are considered in order to choose the appropriate modeling and control approaches to use in the hierarchical control problem.

Components of the DC microgrids

Hereinafter we concentrate on describing the relevant DC microgrid components, namely, their dynamics and the subsequent modeling choices which are made.

Distributed energy sources: Two types of renewable sources are most commonly considered in DC microgrids:

1. Solar panels: Photovoltaics are one of the most commonly used renewable sources in DC microgrids.

PVs generate electricity using solar energy, which is abundant and costless [START_REF] Iovine | Nonlinear control of a dc microgrid for the integration of photovoltaic panels[END_REF]. Their performance depends on the external temperature, the solar irradiation and their location. An important advantage is that they can be connected directly to the DC transmission network through a DC/DC converter;

2. Wind turbines: Wind farms convert the wind energy to electricity through the rotation of the blades [START_REF] Mariam | Microgrid: Architecture, policy and future trends[END_REF]. Wind turbines produce voltage with a variable frequency and AC/DC converters are indispensable to connect with the DC network.

Energy storage devices: As aforementioned, the ES in DC microgrids is very important and necessary for the efficient and continuous operation of the system. The ES system contributes to the cost minimization and the autonomous operation of the microgrid (islanded mode) by storing energy to use when the electricity price is high or when the UG is unable to generate power. Different types of storing devices have been proposed through the years and their development is still in progress. Some of them are mentioned below:

1. Lead acid batteries: This is the most popular type of batteries because they store large amounts of energy and they use a well-established technology which is less costly and more reliable [START_REF] Mutarraf | Energy storage systems for shipboard microgrids-a review[END_REF]. Lead acid batteries are widely available in several sizes and shapes, they are able to handle overcharging situations and they can withstand low and high temperatures. However, some drawbacks exist which influence their capacitance and lifetime such as the number of discharging cycles or the depth of discharges [START_REF] Pavlov | Lead-acid batteries: science and technology[END_REF];

2. Lithium-ion batteries: Recently, the technology of lithium-ion batteries has been ameliorated and they are lately used in DC networks due to their high energy density rate [START_REF] Ma | Control of a hybrid ac/dc microgrid involving energy storage and pulsed loads[END_REF]. Nevertheless, their technology is still under development and changes continuously, a fact which affects their robustness and reliability. Additionally, their intolerance in overcharging or fast discharging and their high cost are significant disadvantages;

3. Ultra-capacitors: Apart from the batteries, ultra-capacitors can be also used since they can deliver current rapidly. Their charging and discharging is fast and they are tolerant of multiple charging and discharging cycles [START_REF] Kanchev | Energy management and operational planning of a microgrid with a pv-based active generator for smart grid applications[END_REF]. Consequently, they can be used for local real-time control to restore the operation of the system in case of emergencies (e.g., blackouts). One major disadvantage is their high self-discharge rate which affects a possible long-term energy storage. Therefore, they are often combined with other types of batteries [START_REF] Morstyn | Cooperative multi-agent control of heterogeneous storage devices distributed in a dc microgrid[END_REF]; 4. Flywheels: This storing device is composed by a motor and uses rotational speed to store kinetic energy.

Although their energy density is high and their lifetime doesn't get affected by the discharging cycles, they suffer from important friction losses [START_REF] Yoldaş | Enhancing smart grid with microgrids: Challenges and opportunities[END_REF].

5.

Vanadium redox batteries: Longevity, fast charging and discharging response, high capacitance, high energy density are some of the advantages of this type of batteries. However, this technology is very recent and the costs of maintenance and installation remain still very high [START_REF] Yoldaş | Enhancing smart grid with microgrids: Challenges and opportunities[END_REF][START_REF] Zubieta | Are microgrids the future of energy?: Dc microgrids from concept to demonstration to deployment[END_REF].

All the aforementioned energy storing devices generate DC current and are connected to the DC microgrid through DC/DC converters for voltage regulation. DC/DC converters: They are electronic regulators which convert the input voltage value to the requested output value. In the industry, two types of regulators exist, linear and switching-based. The linear regulator topology uses linear components (e.g. resistors) to convert the input voltage. Although it is a simpler, less noise-affected electrical circuit, it lacks in efficiency and only the down-conversion operation is possible (higher to lower voltage conversion). On the contrary, the switching regulators use switches to regulate the voltage through a collection of electronic elements like capacitors, inductors and others. Despite their complicated architecture, they are more efficient with lower power losses and have multiple operating modes (Rohm Semiconductors Company, 1997). The main types of switching DC/DC converters are referred below:

1. the buck converter which adjusts a higher voltage to a lower voltage; 2. the boost converter which steps up the input voltage value; 3. the buck-boost converter which implements both of the previously mentioned actions generating a higher or a lower voltage at the output in comparison with the input value.

The majority of the DC/DC converters are unidirectionally supplying voltage from the input to the output. For instance, the Ćuk converter [START_REF] Smedley | Dynamics of one-cycle controlled cuk converters[END_REF]) is a unidirectional buck-boost DC/DC converter which can be connected either to the renewable sources or the loads for voltage adjustment. However, the switching topology allows the possibility of bidirectional functionalities. Such an example is the Split-Pi converter which is a buck-boost bidirectional DC/DC converter and, for instance, can be connected to the ES system for voltage regulation in both directions during charging and discharging [START_REF] Crocker | Power converter and method for power conversion[END_REF][START_REF] Zafeiratou | Dynamical modelling of a dc microgrid using a port-hamiltonian formalism[END_REF]). DC loads: There are two types of loads, the loads which require alternative current (AC loads) and the loads supplied by direct current (DC loads). DC loads have been immensely developed lately due to the rapid evolution of power electronics. Mobile phones, computers, laptops, printers are some of the DC loads used in residential or commercial environments. Besides, LED lighting is another recent technology which operates in DC networks. Finally, the EVs represent a very important chapter of DC loads often considered in DC microgrids [START_REF] Locment | Modeling and simulation of dc microgrids for electric vehicle charging stations[END_REF].

Modeling methodologies

Different modeling methodologies have been already proposed in the literature to describe such complex dynamical systems. The selection of the proper modeling approach depends on the components and the research objectives, i.e. energy management, voltage/current regulation, power balancing, ES sizing and the like. Some modeling methodologies employed through the years are cited below: the Takagi-Sugeno fuzzy modeling approach [START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF] which describes linear models based on input/output datasets covering multiple operation conditions of the system. For instance, in [START_REF] Sáez | Fuzzy prediction interval models for forecasting renewable resources and loads in microgrids[END_REF], the authors use fuzzy modeling to predict the power generation of the renewable source and the loads' demand; the multi-agent based modeling [START_REF] Karavas | A multi-agent decentralized energy management system based on distributed intelligence for the design and control of autonomous polygeneration microgrids[END_REF][START_REF] Khan | Multi-agent based distributed control architecture for microgrid energy management and optimization[END_REF] which concerns agents operating as active entities, individual or collective, and representing a computer model. Using the multi-agent paradigm in a microgrid system, different types of agents can be employed: control agents for controlling the physical units of the system; management agents for managing the microgrid and take decisions; ancillary agents for performing tasks like communication and storage of data [START_REF] Mahmoud | Modeling and control of microgrid: An overview[END_REF]. The global dynamics of the system emerges from the interaction among them; the differential equations modeling [START_REF] Parisio | Stochastic model predictive control for economic/environmental operation management of microgrids: An experimental case study[END_REF], which gives an explicit representation of the system dynamics through constitutive equations. However, through this model description, the power conservation property and the components interconnections are not explicitly deduced; the port-Hamiltonian modeling [START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF][START_REF] Pham | Port-Hamiltonian model and load balancing for DC-microgrid lift systems[END_REF][START_REF] Schiffer | Stability of a class of delayed port-hamiltonian systems with application to microgrids with distributed rotational and electronic generation[END_REF][START_REF] Zafeiratou | Dynamical modelling of a dc microgrid using a port-hamiltonian formalism[END_REF] which generates dynamical models of physical systems taking into account the power-conserving interconnections among the components. The advantage of this method is that the power exchange among the elements is given explicitly and can be used directly for analysis, simulation and control. Besides, PH systems representations can be derived in a straightforward way from their Bond graph models. PH modeling will be employed in this work in combination with Bond graphs and its advantages will be further analyzed in Chapter 2.

Furthermore, the reduction of the model is usually taken into account as a forward step after all the aforementioned methods because of the system's complexity. Sometimes the structure of a system has to be simplified [START_REF] Parisio | An mpc-based energy management system for multiple residential microgrids[END_REF][START_REF] Wang | Intelligent dc microgrid with smart grid communications: Control strategy consideration and design[END_REF] and replaced by simpler dynamical system in order to proceed to its study and supervision. As a next step after the application of the modeling method, the optimal functioning levels of the system must be defined.

Reference profile generation at the upper level control

The reference profile generation aims at determining the optimal operating points in order to respond to the requirements of the system. The optimal profiles are produced providing the optimal power flow over a set of objectives, i.e., cost minimization, power losses reduction, operation of energy storing devices, power balancing. Several methods have been developed for optimal profile generation to deal with the complex dynamics, the nonlinearities, the uncertainties and the constraints. Some of these methods are presented below:

Rule-based controllers take into account a set of statements and limitations which concern the energy management problem of the system. For instance, in [START_REF] Fossati | A method for optimal sizing energy storage systems for microgrids[END_REF] the authors study the capacity of the energy storage system for cost minimization and load forecasting. Furthermore, an energy management algorithm is developed in [START_REF] Moghimi | Rule-based energy management system in an experimental microgrid with the presence of time of use tariffs[END_REF] to generate forecasting profiles for a DC microgrid optimizing the operation of the energy storage system and extracting the maximum power possible from the PV source. The advantage of this method is the easy interpretation which is based on a collection of generic rules under which the system must operate. However, to gather all the possible rules is challenging especially for complex dynamics. Furthermore, the detailed dynamics of the system is not considered and, hence, there are scenarios which cannot be studied.

Mixed integer programming (MIP) is an optimization algorithm which uses the systems' variables as integers. The objective functions and the constraints are presented with linear formulations. In [START_REF] Sechilariu | Power management and optimization for isolated dc microgrid[END_REF], a MIP algorithm is applied to solve the power balancing problem over the simple dynamics of a DC microgrid composed from PV arrays, ES, microturbines and DC loads. The simplicity of the model gives a fast and feasible solution. Noteworthy, the detailed description of the system and the nonlinear constraints are not considered which can lead to inaccurate or unsatisfactory results [START_REF] Braunschweig | Software architectures and tools for computer aided process engineering[END_REF].

Economic Model Predictive Control (MPC) is a closed-loop controller which predicts the behavior of the system under a set of constraints for discrete-time models [START_REF] Rawlings | Model predictive control: Theory and design[END_REF]. Considering a finite horizon, the controller generates at each time slot an optimal open-loop trajectory profile.

Then, only the first step is applied as input to the system to recalculate the states and regenerate the next optimal profile. In [START_REF] Parisio | Stochastic model predictive control for economic/environmental operation management of microgrids: An experimental case study[END_REF]; Prodan and Zio (2014), the authors implement economic MPC for on-line optimal profiles generation to solve the energy management problem.

Differential flatness is an advantageous method used for the analysis of the nonlinear dynamics of the system. To proceed to the optimal profiles generation, it is necessary to isolate concrete physical quantities of the system such as the voltage, the current, the power and others. Therefore, differential flatness can be considered as a suitable tool to inverse the system dynamics, generating in this way all the states and inputs in function of the flat outputs and reducing the number of differential equations. Furthermore, the controllability of the system can be ensured since, according to [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF], a nonlinear system is flat if and only if it is controllable. Recently, differential flatness has been applied in DC microgrids [START_REF] Soheil-Hamedani | Flatness-based control method: A review of its applications to power systems[END_REF]. For instance, in [START_REF] Pham | Port-Hamiltonian model and load balancing for DC-microgrid lift systems[END_REF] the flat representation is applied to generate the optimal profiles for the velocity of an elevator system, while minimizing the energy dissipation. A main advantage of differential flatness is that the resulting profiles implicitly (by their construction from the flat output) respect the system's dynamics. In addition to this, in combination with B-splines approximations, continuous-time constraint validation can be ensured. This method will be further explained and implemented in Chapter 3.

Control approaches

The previous subsection introduced some methods employed in the third level of hierarchical control for optimal power profile generation forecasting the behavior of the system. In general, a control strategy for DC microgrids must be implemented to achieve [START_REF] Unamuno | Hybrid ac/dc microgrids-part ii: Review and classification of control strategies[END_REF]:

Optimization, as previously mentioned, for optimal power flow and cost minimization;

Power quality obtained by a proper power transmission among the sources and the consumers;

Power balancing meaning appropriate load sharing among the components;

Reliability by protecting the microgrid from faulted events or blackouts;

Stability by bounding the error dynamics and ensuring voltage regulation. Therefore, in here, control approaches mostly used at the second or the first level for voltage regulation, power transmission or error mitigation will be presented.

Droop control in DC microgrids is a primary control method which decreases linearly the voltage while the current at the output increases, without ensuring accurate current sharing. For example in [START_REF] Lu | An improved droop control method for dc microgrids based on low bandwidth communication with dc bus voltage restoration and enhanced current sharing accuracy[END_REF], the voltage restoration in the DC bus is achieved over a droop controller implementation. However, the droop control method is not suitable for nonlinear models and its performance in voltage adjustment with accurate load sharing is weak [START_REF] Guerrero | Advanced control architectures for intelligent microgrids-part i: Decentralized and hierarchical control[END_REF].

Plug and play methods (PnP) include information about the subsystems of the microgrid. In this way, the subsystems are able to work independently, meaning that for each component that is added or removed, the stability of the system and the constraints validation are reverified. Furthermore, the controllers of the affected subsystems are readjusted [START_REF] Stoustrup | Plug & play control: Control technology towards new challenges[END_REF]. For instance, in [START_REF] Tucci | Line-independent plug-and-play controllers for voltage stabilization in dc microgrids[END_REF], they aim to prove stability of the linear system using PnP control for voltage regulation. Furthermore, in [START_REF] Riverso | Plug-and-play fault detection and controlreconfiguration for a class of nonlinear large-scale constrained systems[END_REF], the authors use the PnP method in combination with MPC to handle nonlinear networks under constraints and detect faulty behaviors.

Tracking MPC is used for reference tracking. The goal is to alleviate the fluctuations between the reference and the actual profiles satisfying a set of constraints. Tracking MPC is usually combined with economic MPC. Economic MPC generates the optimal reference profiles which, subsequently, are sent to the tracking MPC as inputs to follow them under perturbations [START_REF] Legry | A two-layer model predictive control based secondary control with economic performance tracking for islanded microgrids[END_REF]Zafeiratou et al., 2018a). Tracking MPC will be implemented in this work at the second level of the hierarchical control. Additionally, there is the robust MPC variant which builds upon an optimization problem influenced by uncertainties and tends to increase the robustness of the system [START_REF] Dragičević | Model predictive control of power converters for robust and fast operation of ac microgrids[END_REF].

Agent-based controllers exchange information among local entities with limited knowledge. The agents can interact to each other to solve an optimization problem for load balancing, voltage regulation or battery charging. This method provides autonomy and many controllers can be implemented in an agent-based framework. For instance, the authors of [START_REF] Li | Agent-based decentralized control method for islanded microgrids[END_REF] propose a two-layer decentralized control, where the lower layer controls the power distribution while at the upper layer the communication among the agents is applied.

In general, the control methods are divided into two categories, centralized and decentralized. In a centralized supervision, a central controller performs the requested actions to achieve proper power distribution and secure operation of the microgrid. Although, it is an advantageous method because it facilitates the controllability and the observability of the system, the controller can collapse under a single failure. Furthermore, the complexity in computation does not allow its implementation in large-sized microgrids. On the other hand, the decentralized control does not demand global information of the system, since this method is based on the local supervision of a subsystem or a unit. The decentralized controllers do not require real-time communication, but the observability on the global behavior of the system is limited [START_REF] Meng | Review on control of dc microgrids and multiple microgrid clusters[END_REF].

Thesis orientation

The aforementioned literature demonstrates a wide variety of modeling and control approaches performed in DC microgrids. This thesis will concentrate on the analysis of a meshed DC microrgid proposed in the ANR C 3 µ project (Components, Control, Communication) (Agence Nationale de la Recherche, 2016). The project aims at optimizing the power distribution in commercial or residential buildings employing a central management system. The goal is, firstly, to satisfy the consumers' needs associated with an effective power generation from the renewable energy sources and the ES system. The originality of the project resides upon the definition of a smart node (SN), which is composed by two Split-Pi DC/DC converters explicitly described in Chapter 2. The SN will manage the central transmission network and the power exchange among the components as well as the reconfiguration of the system. The main tasks of the project are summarized below:

Design and implementation of the DC/DC converters and the components;

Modeling of the EMS and study of the reliability of the DC microgrid; Verification under simulations to estimate the proposed method's viability for the power balancing problem;

Design of the communication network among the SNs.

The end goal is to remove the central EMS and to integrate all the implemented actions in the SNs, which will supervise directly the DC microgrid's operation.

In this manuscript, the energy management problem is considered. The main objective is to formulate a hierarchical control problem which will optimize the power distribution within the DC microgrid. It will also provide a solution for the voltage regulation and the discrepancies' mitigation among the optimal and the real profiles. The power flow will be optimized regarding the electricity cost, the consumers' demand, the power dissipation and the ES system. To handle all these aspects, the following must be considered: the switching control of the converters is necessary. To achieve this, the detailed knowledge of their physical characteristics (capacitors, inductors, resistors and so on) is indispensable; the number of converters must be limited to prevent an increase in power losses. Note that the converters will handle both the power flow and the voltage/current regulation of the system; the best possible exploitation of the renewable energy sources is important. The unused energy must be either stored or sold to the main grid; a collection of batteries must be considered which will ensure the safe and reliable operation of the system in case of faults or uncertainties; a specific amount of transmission lines must be placed in a way that the power will be distributed under the less potential power losses. The lines under fault could be also isolated. Therefore, it is clear that the overall problem has many aspects and many problems remain open. The main target is to have a unique dynamical model where every physical component will be explicitly described. Port-Hamiltonian (PH) representations give this flexibility of describing the dynamical model in detail. The DC microgrid system will be, firstly, presented in PH form giving high importance to the ES system. Next, the computation of the flat outputs of each component will be provided, generated by an already-proposed algorithm in [START_REF] Franke | On the computation of flat outputs for nonlinear control systems[END_REF]. In addition, the flat representations of the PH systems will be further analyzed under the bicausality notion which investigates the inverse dynamics of the system. Afterwards, the constrained optimization-based control approach will be introduced by employing the differential flatness method in combination with a B-splines parametrization.

In the sequel, the hierarchical control approach will be implemented as a three-levels scheme: the high, the middle and the low level. At the high level, optimal profiles will be produced to be followed by the middle and lower levels, respectively. Power balancing requirements at the high level will lead to an optimal profile generation for the battery, to be tracked at the middle level. Then, the middle level will provide voltage and current references for the battery. These profiles will be tracked at the low level by an explicit switching law used in the DC/DC converters. The different timescales, the nonlinearities and the constraints of the system will be also included firstly for the cost minimization. Afterwards, the power losses of the central transmission network will be considered to generate an optimal power flow. The reference profiles generated by differential flatness and the B-splines parametrization will be compared with the optimal profiles obtained by economic MPC. Furthermore, other scenarios will be taken into account as for example for different distances among the power sources or when lines under fault appear.

Contributions of the thesis

This thesis presents a meshed DC microgrid architecture in PH representation supervised by a hierarchical control strategy to manage load balancing and efficient power distribution. A novel combination of PH modeling with differential flatness and B-splines parametrization is introduced and leads to an off-line constrained optimization problem. This method provides continuous-time optimal profiles not only for the power, but also for the voltage, the current and other physical quantities of the electrical circuits (the charge, the magnetic flux and the like). Hence, it can be an advantageous and reliable solution to study thoroughly and predict the behavior of the system.

More precisely, the main contributions of this work are summarized below:

the dynamical model representation with PH formalism through its associated Bond graph: PH method generates a well-structured dynamical representation of the system. This method naturally provides power-preserving differential equations. Additionally, its structure allows the isolation of any physical quantity for investigation (voltage, current, charge, magnetic flux) without changing the primary mathematical model. The proposed meshed DC microgrid is presented in PH form. Firstly, the PH state-space representation is given for every part of the microgrid separately. Next, the PH state-space representation of the central transmission network is analytically introduced; the algorithm proposed in Franke and Robenack (2013): the algorithm is implemented for the PH system to find its flat representation. Furthermore, the bicausality notion is applied to further analyze the flat representation of the PH models and study a more direct way to calculate the set of flat outputs. Comparisons amongst the different set of flat outputs are realized through simulations. The goal is to find a straightforward way to calculate the flat outputs of the PH systems. A method is proposed for electrical circuits with illustrative examples; the use of differential flatness and B-splines parametrization for optimal profile generation: microgrids are complex systems with complicated mathematical models and numerous states. In this work, differential flatness seems to be an appropriate tool which can inverse the dynamics of the system and describe the states in function of the flat outputs. Furthermore, using B-splines to parametrize the flat output allows, due to the properties of B-splines, to characterize the system's functioning in continuous time. Optimal reference profiles for the UG and ES power will be obtained. In addition to the previously referenced publications, which concentrate more in obtaining optimal profiles for the power or the energy generated from each source, in here we go further and produce optimal profiles in continuous-time not only for the power but also for the current and the voltage of the components;

the proposed hierarchical control approach to solve the power balancing problem: works in the literature concentrated mostly on one control level. In the case where they used more than one level, they simplified either the model or the control method. In this work, a complete hierarchical control framework is presented which investigates many important factors at the same time: cost minimization, power balancing, power losses, constraints, tracking errors (decreasing the voltage and current discrepancies) and regulation of the switches in converters.

We provide, hereinafter, the list of publications accepted/submitted to various conferences and journals, and the presentations given by the candidate at various scientific days or mobility projects:

Organization of the manuscript

This thesis consists of 5 chapters (see Fig. 1.4.1 on page 11), including the introduction.

Chapter 2, firstly, presents some notions and definitions for the Bond graphs and the PH systems. Next, the Bond graph of each component of the meshed DC microgrid is developed. The associated PH state-space representations are provided. Furthermore, a simplified version of the central transmission network is presented. The transformation from the complex to the simplified version is described through their PH formulations. Finally, the constraints and the control objectives of the system are introduced;

Chapter 3 presents the flat representation of each component of the system. The detailed explanation of the algorithm's implementation is provided taking as an example the ES PH representation. Simulations for the verification of the computed flat outputs are also included. Then, the theory of bicausality is introduced. The flat representation through bicausality is presented through illustrative examples of some electrical circuits of the system; Chapter 4 describes the hierarchical control problem. Firstly, some basic tools of the hierarchical control problem formulation are provided. Afterwards, the energy management problem is presented together with the control objectives, the reference profiles and the set of constraints. Next, the multilayer supervision problem is analyzed explicitly with extensive simulation results which are explained in detail. Furthermore, a comparison with MPC as another method of optimal profile generation is implemented. Finally, different scenarios for the DC microgrid are investigated (e.g. lines under fault); Chapter 5 draws the conclusions of the thesis and analyzes the short-term and long-term perspectives.

DC microgrid dynamical representation 2.1 Introduction

Through the years, microgrids have become progressively convoluted power systems since they are composed by a large number of small-scale subsystems such as a collection of renewable energy resources, an enormous ensemble of power electronics converters, a set of numerous loads, an assortment of batteries for energy storage (ES) and a multi-line transmission network. Therefore, an appropriate structure plays a crucial role for an effective operation and, towards this direction, there are many aspects to investigate. For this reason, results from the state of the art implement mathematical models to represent and study the dynamics of such systems. The complexity of these networks requires a concrete modeling methodology according to the objectives, some of which are enlisted below: which microgrid type to use, AC, DC or a combination of both, depending on the distributed energy resources (DERs) included, the kind of loads, the amount and the category of converters; the selection of the transmission lines and the network topology, radial or meshed, the number of transmission lines and the distance among the DERs and the loads. These characteristics rely upon the financial planning of the consumers, the cost of the energy purchase from the main grid, the power dissipation and the like; the energy management method, centralized, decentralized or distributed. Each approach needs a different analysis, meaning a different modeling method and control strategy; the analysis of the power flow, signifying, primarily, the in-depth study of the variables of the electronic devices included in the system, such as the input/output voltage of the sources, the switches of the power converters, the analysis of the fast and the slow scale dynamics of each part; the necessity of the ES, which influences the reliability and the efficiency of the system. The sizing of the batteries, the depth of discharge, the deterioration and the lifetime are some of the factors that influence the cost and the performance in such systems; which renewable energy sources to use, depending on their type, AC or DC sources, current or voltage sources. The buying cost must be also considered not only for the sources themselves but also for the converters needed for their synchronization with the central transmission network. This work, as already mentioned, is part of the ANR-15-CE05-004-02 C 3 µ project (Components, Control and Communication) which proposes a meshed DC microgrid architecture to deliver electrical energy to the interior of a building (commercial or residential) under a decentralized energy management. The project is investigating the energy saving potentials, the benefits, the barriers and the efficiency of employing directly the DC current for local power distribution rather than the AC current which needs to be converted for the supply of the DC loads. The main target is to have a well-structured, explicit dynamical model which respects and follows the constitutive laws of physics and engineering. Port-Hamiltonian (PH) representation offers this flexibility by describing a physical system in detail and by giving information separately for each 13 physical quantity, such as the charge, the voltage, the dissipation and the like, to use it afterwards for further analysis.

In the literature, PH modeling has been considered for various dynamical systems. For instance, in 2014)). PH modeling has been proposed even for economical modeling in [START_REF] Machado | A fractional perspective to the bond graph modelling of world economies[END_REF]. Moreover, in [START_REF] Macchelli | Port-hamiltonian formulation of simple macro-economic systems[END_REF], the authors take advantage of the energy storing and dissipative elements, the power exchange and power conservation in order to analyze the market's supply and demand.

In general, PH modeling is a structural approach which focuses on providing an exhaustive, explicit and modular description of the power routing in a system through a complex port-based network. The PH statespace representation represents the constitutive equations of its energy storing elements coming from various domains, chemical, mechanical, electrical, thermodynamical and the like. Furthermore, it can be derived in a straightforward way from its associated Bond graph model, a graph-oriented approach which leads to the direct acquisition of the PH state-space representation, as will be explained later in the text [START_REF] Karnopp | System dynamics: modeling, simulation, and control of mechatronic systems[END_REF].

In this chapter, the meshed DC microgrid system will be presented composed by: an external utility grid (UG), a set of solar panels (PV), a series of electrical vehicles (EV), a collection of lead-acid batteries and a group of loads. DC/DC converters are also included for voltage regulation and for power flow direction. The main idea of this chapter stems from the detailed representation of the DC microgrid describing each component separately, presenting its Bond graph and its associated PH formulation.

In the sections that follow, the main components of the DC microgrid will be presented. Afterwards, we will introduce the basic theory and definitions of the Bond graph representation and the PH formalism. Then, we will apply these methods in the DC microgrid replacing the components (sources, converters, transmission network) with electrical circuits. Next, after introducing the complete version of the DC microgrid, we will develop a reduced version preserving the primary characteristics. In this version, the ES system will play an important role and will be described explicitly. At the end, the contributions of this chapter will be mentioned. Note that, for the sake of convenience, wherever it is straightforward implied by the text, the time dependence will be ignored.

DC microgrid architecture

The meshed DC microgrid considered in the project (Fig. 2.2.1), as already mentioned, is composed by a renewable energy source (PV), an ES system composed from lead-acid batteries and a number of loads (office and housing devices, electrical vehicles (EV) and the like). The microgrid is connected to a three-phase UG through a bi-directional AC/DC inverter and with the PV, they constitute the main power sources of the system. Between the UG and the DC microgrid, there is a DC breaker with a switch which disconnects the DC microgrid from the AC-grid avoiding occasional blackouts or other fault events. When the DC breaker's switch is open, the DC microgrid is able to operate autonomously in islanded mode. The common DC bus is a high-voltage transmission network of 400 V DC.

In the following, the switching DC/DC converters of the system will be presented. These components are regulating the input/output variables among the sources, the loads and the DC bus. Furthermore, they control the power flow direction within the central transmission network. In this work, two types of converters will be used: a bi-directional buck-boost converter, called Split-Pi converter [START_REF] Crocker | Power converter and method for power conversion[END_REF], which produces an output voltage higher or lower than the input voltage. At the same time, this converter can operate in both directions which makes it suitable for sources which generate and accept energy, such as the battery (charging/discharging). To implement these parallel functionalities, four switches exist, which are the main components responsible for its efficient operation. A pair of Split-Pi converters constitutes the smart node (SN) (see also Fig. 2.2.1) of the system which takes decisions related to the power demand, the electricity cost, the load balancing and the state of charge of the ES. In other words, this an unidirectional buck-boost converter, called Ćuk converter [START_REF] Pyrkin | Ćuk converter full state adaptive observer design[END_REF]. This converter regulates the output voltage which can be lesser or greater than the input. It is appropriate for sources or loads which either generate or receive power.

A significant factor of the converters is the duty cycle1 , d, of the switches which characterizes the converters' activity. In the following, the aforementioned components of the microgrid will be described. As mentioned before, the main renewable source of the system is the solar panel. The solar panel is composed by arrays of solar cells in parallel or in series. The solar cell is an electrical device that collects the solar irradiation and generates electricity. A classical solar cell can be illustrated with a current source and a diode in parallel. Since the solar cell is not an ideal element without power losses, two resistances are added in the circuit, a shunt resistance R 1pv and a series resistance R 2pv [START_REF] Villalva | Comprehensive approach to modeling and simulation of photovoltaic arrays[END_REF]. The shunt resistance decreases the amount of the current passing through the electrical circuit. On the other hand, the series resistance lowers the quality of the solar cell by reducing the fill factor2 . In Fig. 2.2.2, the equivalent electrical circuit is depicted. The PV panels produces direct current and they are appropriate renewable sources for a DC network. In the literature, several models have been proposed which describe the battery's output voltage and the operations during charging, discharging and recovery. The model selected for this work is based on the chemical kinetics process and is called Kinetic Battery Model (KiBaM) [START_REF] Manwell | Lead acid battery storage model for hybrid energy systems[END_REF]. It is considered as a two-tank model, the bound charge and the available charge, separated by a conductance k b (Fig. 2.2.3). The tank q 1b supplies electrons directly to the load and the tank q 2b supplies electrons only to the available charge tank [START_REF] Jongerden | Battery aging, battery charging and the kinetic battery model: A first exploration[END_REF].

Description of the meshed DC microgrid components

Lead acid battery

The flow between the two tanks depends on their height difference and on the parameter k b , which is constant. As we can see from fig. 2.2.3, the two heights h 1b and h 2b are given by the following equations (Jongerden and Haverkort, 2017):

h 1b (t) = q 1b (t) c b , (2.2.1)
h 2b (t) = q 2b (t) 1 -c b (2.2.2)
with the charging variations equal to:

dy 1b dt = -i b (t) + k h 2b (t) -h 1b (t) , (2.2.3
)

dy 2b dt = -k b h 2b (t) -h 1b (t) , (2.2.4)
with initial conditions y 1b = c b C b and y 2b = (1c b )C b , where q 1b is the available charge, q 2b is the bound charge, c b is the width of the tank q 1b , which refers to the total capacity of the available charge, 1c b is the width of the tank q 2b , which refers to the total capacity of the bound charge, and C b is the total capacity of the battery. Moreover, i b is the current of the battery delivered to the load. When a load is applied to the battery, the available charge decreases. When, the load stops consuming, the bound charge starts to flow in the available charge tank until the heights are equal again.

In the following, the KiBaM model is replaced by an equivalent electrical circuit [START_REF] Pham | Port-Hamiltonian model and load balancing for DC-microgrid lift systems[END_REF]. Considering the two tanks as two storage elements, C 1b and C 2b , connected together with a resistance, R 2b , we obtain the electrical circuit depicted in Fig. 2.2.4. Furthermore, a resistor, R 1b , is added to represent the transmission line between the source, v sc , and the battery.

Split-Pi converter

Next, the Split-Pi converter and its operation are presented. This converter is bi-directional and produces an output voltage which can be higher or lower than the input voltage in both directions. It is a buck-boost converter (DC/DC converter) and it is controlled by four switches. It has four modes of operation: two up-conversions (higher output voltage than the input) and two down-conversions (lower output voltage than the input) [START_REF] Crocker | Power converter and method for power conversion[END_REF][START_REF] Singhai | Modeling and analysis of split-pi converter using state space averaging technique[END_REF]. We start by presenting in Fig. 2.2.5 the electrical circuit of the Split-Pi converter and in Fig. 2.2.6 its operation modes as presented also in Table 2.2.1. We mention at this point that we consider the switches as ideal, lossless elements, which can conduct current at zero voltage, when Sw is equal to 1 and hold a voltage at zero current, when Sw is equal to 0 (Escobar et al., 2015; van der Schaft and Jeltsema, 2014). Hereinafter, we continue with the Ćuk converter which is unidirectional, based on a buck-boost topology and produces voltages which can be higher or lower than the input voltage [START_REF] Escobar | A hamiltonian viewpoint in the modeling of switching power converters[END_REF]. This kind of converter regulates the input voltage of the loads and the EVs and the output voltage of the PV system.

In the following figures, we present the electrical circuit of the Ćuk converter (Fig. 2.2.7) and its modes of operation (Fig. 2.2.8).

In the sequel, the mathematical model of the aforementioned components will be presented. The Bond graph method will be employed in order to develop the PH state-space representations of the system. The Bond graph is a graphical representation which describes the energy-preserving structure of a dynamical system. This method, built upon the power interconnections amongst the physical elements, clarifies the system's behavior and contributes to the formulation of the system's equations.

Primarily, the Bond graph is regarded as an interconnection of three types of components, as you can see also in Fig. 2 With these three components, Bond graphs cover multiple domains and describe various systems, such as mechanical, electrical, thermodynamical and others. Regarding the multi-scale systems, using the notion of Bond graph, they can be decomposed and analyzed as reduced interconnected parts. The strokes in the Bond graph represent the power flow through the elements with a half-arrow pointing the power flow direction. Each arrow is characterized by a pair of effort, e, and flow, f , and their product is the power passing through, P = e • f . The components are connected together in the Junction Structure, JS (Fig. 2.3.1). In Table 2.3.1, several interpretations of e and f are depicted for different types of systems. Afterwards, the main characteristics [START_REF] Borutzky | Bond graph based physical systems modelling[END_REF] of the Bond graph will be explained through an example of an electrical circuit (Fig. 2.3.2): according to Fig. 2.3.2, the electrical circuit has two storing elements, one capacitor C and one inductor I, one dissipative element R and two sources, one voltage source S e and one current source S f . These elements are considered as one-port elements referring to one single port and one pair of e and f ; in Fig. 2.3.3 the corresponding Bond graph of the electrical circuit is presented. The arrows demonstrate the power flow, P = e • f , within the system and compose the Junction Structure; in the Bond graph 0 and 1 junctions also exist. Their major characteristics is the power conservation and the reversibility. 0 junction has equal efforts and the sum of flows is equal to 0 (e 1 = e 2 = e 3 = ... = e n , f 1 + f 2 + f 3 + ... + f n = 0). On the other hand, 1 junction has equal flows and the addition of flows sums up to 0 (f 1 = f 2 = f 3 = ... = f n , e 1 + e 2 + e 3 + ... + e n = 0). This is equivalent to the rules which the electrical circuits follow in parallel and in series respectively; finally, the storing and dissipative elements obey to the following constitutive equations: the R element

respects e R = f R • R and P = e 2 R R = f 2 R • R, C follows e C =
q C C and I obeys to f I = p I I . In addition to the previous statements, the 0 and 1 junctions follow also specific causality rules. In 0 junctions the effort is established and there is only one causal stroke with the causality symbol pointing towards the junction as in Fig. 2.3.7(a). Whereas, for 1 junctions, the flow is established, hence there can be only one causal stroke with the causality symbol pointing towards the element as in Fig. 2.3.7(b). For better comprehension, the causal Bond graph of the aforementioned electrical circuit(Fig. 2.3.2) is provided as well as the procedure of the causality assignment. Furthermore, two more variables must be explained, which are necessary for the formulation of the system's state-space representation. These are the momentum, p, and the displacement, q, variables. The p variable is considered to be the e time integral (p(t) = e(t)dt) and the q variable is the f time integral (q(t) = f (t)dt). In the case of electrical circuits, p represents the magnetic flux of the inductors and q the charge of the capacitors. Therefore, according to the integral causality, we further obtain:

f or C : e = q C , f or I : f = p I .
Therefore, from the causal Bond graph of the electrical circuit example (Fig. 2.3.10), e C and f I can be replaced by:

f or C : e C = q C C , f or I : f I = p I I .
In the aforementioned description, only the elements which will be used in this work were presented. Although, in a Bond graph, two-port elements can also be considered (transformers, T F , and gyrators GY ). These elements do not create or store energy. They transmit power by changing the e f ratio (Fig. 2.3.11).

Port-Hamiltonian systems

In this section, we introduce the port-Hamiltonian (PH) systems together with some basic definitions. Additionally, the conversion from Bond graphs to port-Hamiltonian state-space representations will be presented. These notions will be instrumental for the mathematical model representation of the DC microgrid of Fig.

ẋ(t) = [J(t) -R]Qx(t) + Gu(t), y(t) = G Qx(t) + Du(t). (2.3.1)
This is a typical PH representation and each of the components has a clear physical justification:

x(t) ∈ R n is the state vector of the extensive energy variables p(t) and q(t), referred before as the displacement variables and derived from the energy storing elements; u(t) ∈ R m is the input vector, which contains voltage or current variables, generated by the sources or received by the loads; y(t) ∈ R m is the output vector, which contains voltage or current variables, generated by the PH system; J(t) ∈ R n×n describes the interconnections of the system and defines the power continuity, meaning the structural balance equations which typically contain no numerical parameters. It is a skew-symmetric and positive semi-definite matrix;

R ∈ R n×n is likewise a skew-symmetric positive semi-definite matrix which outlines the dissipative elements;

Q ∈ R n×n is a positive definite matrix which, in the linear case, is a diagonal matrix including the system's physical parameters, such as the capacitance C of the capacitors and the inductance I of the inductors. More complex non linear constitutive equations may also be included;

D ∈ R m×m describes the direct interconnection of the input variables;

G ∈ R n×m is the control matrix of the PH system.

The output vector y(t) in (2.3.1) is determined according to the inputs selection and the input map G. The inner product among the input and the output vectors minus the dissipative energy gives always the external power supplied to the system [START_REF] Escobar | A hamiltonian viewpoint in the modeling of switching power converters[END_REF]. Below, we examine the Hamiltonian in linear cases:

1. The Hamiltonian, H, is the total energy stored in the storage elements of the system equal to:

H(x) = 1 2 x(t) Qx(t); (2.3.2)
2. From (2.3.1) and (2.3.2), when D = 0, we obtain the external power supplied to the system:

d dt H = u(t) y(t) -x(t) Q RQx(t) -u(t) Du(t). (2.3.3)
In general, PH systems are determined by the power-preserving interconnections of the elements of a physical system through a Dirac structure (DS) (van der Schaft et al., 2014). As in the Bond graphs with the Junction structures, a DS classifies equally power-conserving interconnections into three types of elements: i) energy-routing elements; ii) energy-storing elements ; iii) energy-dissipative elements. Energy-routing elements are the transformers and gyrators. Additionally, energy-dissipative elements are the resistors and energy-storing elements are the capacitors and the inductors. Besides, in a DS, we have external ports which determine the interaction of the system with its external environment. In Fig. 2.3.12, we observe the correspondence between Bond graphs and DS which both account for the continuous power balance equations (e.g., Kirchoff's laws in electrical networks). The DS connects all the port variables together and obeys the power conservation law which means that the total power e f of the central network is equal to 0 [START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF]van der Schaft et al., 2014). To define a DS, we consider a finite-dimensional linear space of flows F with f elements, called flow vectors. The efforts' linear space will be denoted as E = F * and contains the e variables. The total space is equivalent to F × E and it is called the space of port variables and the total power is their duality product (van der Schaft et al., 2014):

< e|f >= e f, (f, e) ∈ F × E, (2.3.4)
e f = 0, dimD=dimF.
Looking into the geometrical structure of the total space F × E, according to the power definition, the canonically bilinear form of a DS can be defined as follows (van der Schaft et al., 2014):

<< (f a , e a ), (f b , e b ) >>:=< e a |f b > + < f a |e b >, with (f a , e a ), (f b , e b ) ∈ F × E. Definition 2.3.2. (Dirac structure -van der Schaft et al. (2014)): A Dirac structure on F × E is a subspace D ⊂ F × E such that D = D ⊥ ,
where ⊥ denotes the orthogonal supplement with respect to the bilinear form , .

From the DS port-based graph, we can derive the state-space representation of a dynamical system in PH form. The PH system results from the linkage of the DS with the elements associated to the storage, the dissipation and the external environment. 2014)): Consider the total storing energy as the Hamiltonian H defined by a state space X , where H : X R, its tangent space T x X and co-tangent space T *

x X . Then, the PH system is a DS defined by:

D ⊂ T x X × T * x X × F R × E R × F E × E E ,
with the ES ports regarded as (f S , e S ) ∈ T x X × T * x X, for which f S =ẋ and e S = ∇H(x), the dissipative port variables regarded as:

R R = {(f R , e R ) ∈ F R × E R |R(f R , e R ) = 0, < e R |f R >≤ 0}
and the external ports as (f E , e E ) ∈ F E × E E . In general, a PH system is expressed as noted below:

(-ẋ, ∇H(x), f R , e R , f E , e E ) ∈ D.
In the following, the Bond graph theory in combination with the PH state-space representations will be employed to describe explicitly the dynamical model of each of the DC microgrid component mentioned before in Section 2.2.1.

DC microgrid dynamical model

In this section, we will present the complete model of the meshed DC microgrid and introduce its global Bond graph representation. Hereinafter, we will provide the mathematical models of the energy sources, the PV panel and the lead-acid batteries, in PH form based on their Bond graphs. Afterwards, we will describe also the central transmission network and its characteristics. In general, it is composed by three meshes, the PV station with its ES system, the EV station and the various loads (printers, computers, LED lighting, mobile phones and the like). Furthermore, the transmission lines have been replaced with an 'inductor plus resistance' block. This is the so-called RL circuit in which the two components are placed in series. In DC networks, an inductor behaves as a static resistance and contributes to the reduction of the power losses.

Dynamical model of the sources

Before beginning to elaborate separately each part of the DC microgrid, we mention, at this point, that the DC microgrid will be considered and analyzed as an ensemble of electrical circuits, where the storage elements are the inductors, I, and the capacitors, C, the dissipative elements are the resistors, R, and the energy sources are the voltage sources, as has been described in Section 2.3.1. The energy stored in the electrical circuits (Fig. 2.3.2 for instance), otherwise called the Hamiltonian, is a function of the magnetic flux of the inductors, p, and the charge of the capacitors, q. Therefore, from the general description of the Hamiltonian in linear cases (2.3.1), a particularization for the case at hand is given below:

H(x) = 1 2 x(t) Qx(t),
where

x(t) = p 1 p 2 ... p n I q 1 q 2 ... q n C ∈ R n×1 is the state vector (n = n C + n I ,
where n I is the number of inductors and n C the number of the capacitors). Q is the circuit parameter matrix equal to

Q = diag 1 I 1 , 1 I 2 , ..., 1 I n I , 1 C 1 , 1 C 2 , ..., 1 C n C ∈ R n×n .
Therefore, the Hamiltonian becomes:

H(x) = 1 2 p 2 1 I 1 + p 2 2 I 2 + ... + p 2 n I I n I + 1 2 q 2 1 C 1 + q 2 2 C 2 + ... + q 2 n C C n C . ( 2
v s i ph + v Di i Di + v R1pv i R1pv + v R2pv i R2pv = 0, (2.4.2)
where s is the sunlight source, the series resistance R 2pv resists to the current flow, R 1pv is the shunt resistance of the system, Di is the diode and v declines the output voltage. Additionally, from the 0/1 junctions, the equations below are deduced:

i ph = i Di + i R1pv + i 2pv , v R2pv = v Di . (2.4.3)
However, the model of a PV panel is rather complicated and obeys a set of constitutive equations that have been established the past decades from the researchers and can be found in the literature [START_REF] Enrique | Theoretical assessment of the maximum power point tracking efficiency of photovoltaic facilities with different converter topologies[END_REF][START_REF] Villalva | Comprehensive approach to modeling and simulation of photovoltaic arrays[END_REF]. Therefore, all the unknown variables of (2.4.2) and (2.4.3) are calculated by the classical PV mathematical model. The parameters of the PV panel depend on the selected model which, in our case, is the DS-100 MPV module (180 W peak PV generation). The final equation of the PV panel composed by cell arrays is the following:

i(t) = N pa i ph (t) -N pa i 0 (t)     exp     v(t) N se + i R 2pv N pa nv Di (t)     -1     -i 1pv (t), where (2.4.4) v Di (t) = k B T emp q el , (2.4.5) i 1pv (t) = v(t) N pa N se + i(t)R 2pv R 1pv , (2.4.6)
P (t) = v(t)i(t), (2.4.7)
where i is the output current, N pa is the number of cells in parallel, N se is the number of cells in series, i 0 is the saturation current of the cell, k B is the Boltzmann's constant, q el is the electron charge and T emp is the working temperature of the solar cell.

Voltage (V) To proceed to the PV profiles generation, the model was developed in MATLAB/Simulink. For model testing, the aforementioned PV model is studied under a constant input solar irradiation equal to 800W/m 2 and a constant external temperature equal to 25 • C to obtain the I-V curve in Fig. 2.4.3a and the P-V curve in Fig. 2.4.3b.

For further validation, we generate also power profiles (Fig. 2.4.5) based on different weather scenarios (Fig. 2.4.4a) for a PV array of 28 panels with 72 solar cells. Besides, we take into account that a solar cell produces about 0, 5V . A typical PV that produces about 12V contains 36 cells wired in series. For a PV with 72 cells, we consider two 12V PV panels of 36 cells wired in series, usually with a jumper, allowing an output voltage equal to 24V . Therefore, we proceed to the following simulations and the obtained profiles will be used later in Chapter 4. Consequently,from the Bond graph, we arrive directly to the following PH state-space representation:

Time(h)

           q1b q2b = [J b -R b ]Q b x b +   - 1 R 1b 0   (-v sc ), i sc = - 1 R 1b 0 Q b x b + 1 R 1b v sc , (2.4.8)
where x b = q 1b q 2b ∈ R 2×1 is the state vector, v sc ∈ R is the system's input vector denoting the input voltage coming from the Split-Pi converter, presented afterwards, i sc ∈ R is the output vector and where i sc is the battery's current during charging. The dissipation matrix R b ∈ R 2×2 is equal to:

R b =    1 R 1b + 1 R 2b - 1 R 2b - 1 R 2b 1 R 2b    .
(2.4.9)

All the unknown variables and parameters of the battery model can be found similarly in the literature [START_REF] Manwell | Lead acid battery storage model for hybrid energy systems[END_REF] and depend on the type of the lead-acid battery. The charging and discharging processes of the battery will be controlled by the Split-Pi converter, presented in the next section. As it can be deduced by (2.4.8), the matrix J b of the battery PH representation is equal to 0. The matrix J includes the internal power-preserving interconnections representing the oscillations between the electrical and the magnetic field of the storing elements. Since, there is only capacitors in the electrical circuit and no inductors, no oscillations exist and the J b matrix becomes 0. Therefore, by connecting the battery circuit to the Split-Pi converter, as it will be shown later, we obtain the J es matrix of the ES system. For the simulations, an AGM 12-165 battery (165 Ah battery capacity) is considered and is validated later within the ES system (subsection 2.5.1).

Dynamical model of the converters

Split Pi converter model

This subsection introduces the corresponding Bond graph (Fig. 2.4.7) of the Split-Pi converter (Fig. 2.2.5). In the Bond graph, we add a resistor element after the source, R 1sc , in order to respect the causality property (integral causality) and avoid the differential causality (see also the causality discussion carried out in Section 2.3.1). According to the Bond graph of the converter, its PH state-space representation is provided below:

                                     ṗ1sc ṗ2sc q1sc q2sc q3sc         = [J sc -R sc ]Q sc x sc +          0 0 0 0 -1 R 1sc 0 0 0 0 1          u sc , y sc =   0 0 -1 R 1sc 0 0 0 0 0 0 1   Q sc x sc +   - 1 R 1sc 0 0 0   -v s -i R Lsc ,
(2.4.10)

where x sc = p 1sc p 2sc q 1sc q 2sc q 3sc ∈ R 5×1 is the state vector,

u sc = -v s -i R Lsc ∈ R 2×1
is the system's input represented by the input voltage v S generated from the source, the PV or the ES, and the current of the load i matrix R sc ∈ R 5×5 are described below:

R Lsc = q 3sc R Lsc C 3sc , y sc = i s v R Lsc ∈ R 2×1 is the output vector, i S is
J sc =       0 0 1 -(1 -d 1sc (t)) 0 0 0 0 (1 -d 2sc (t)) -1 -1 0 0 0 0 1 -d 1sc (t) -(1 -d 2sc (t)) 0 0 0 0 1 0 0 0       , (2.4.11) R sc = diag(0, 0, 1 R 1sc , 0, 0), (2.4.12)
where values d 1sc (t), d 2sc (t) are the control variables representing the duty cycles of the converter's switches. As a next step, the relations among the duty cycles needs to be defined. This will help us to develop an open-loop control law for the switches. According to the patent of the United States Patent and Trademark Office No: US 6914420 B2 published on July 2005 [START_REF] Crocker | Power converter and method for power conversion[END_REF], the following statements can be deduced for the down and the up-conversion (see also Table 2

.2.1):

During the down-conversion (phase 1), the C 2sc voltage equals the input voltage v S , while the current rate of the inductor

I 2sc is equal to di I2sc dt = v I2sc I 2sc .
During the down-conversion (phase 2), when the voltage around I 2sc is reversed, the current rate of the

inductor I 2sc is equal to di I2sc dt = -v C3sc I 2sc .
Therefore, the absolute values of the two current ratios over a cycle must be equal and for each ratio we can define a period of time, t 4 for the first ratio and t 3 for the second:

t 4 v C2sc -v C3sc I 2sc + t 3 -v C3sc I 2sc = 0 (2.4.13) v C3sc v C2sc = t 4 t 3 + t 4 = α, (2.4.14)
where α is a factor that will help us to determine the switching action within the converter. Now, if we take into consideration that, by definition,

d 2sc = t 3 t 3 + t 4
, then a = 1d 2sc . For the up-conversion, we have similarly:

During the up-conversion (phase 1), the C 2sc voltage equals to the output voltage of the circuit, while the current rate of the inductor

I 1sc being equal to di I1sc dt = v I1sc I 1sc .
During the up-conversion (phase 2), the current rate of the inductor

I 1sc is equal to di I1sc dt = v C1sc I 1sc .
Consequently, as before, the absolute value of the two current ratios over a cycle must be equivalent and for each ratio we define another two periods of time, t 2 and t 1 :

t 2 v C1sc -v C2sc I 1sc + t 1 v C1sc I 1sc = 0 (2.4.15) v C2sc v C1sc = t 1 + t 2 t 2 = α. (2.4.16) If we account that d 1sc = t 1 t 1 + t 2 , then a = 1 1 -d 1sc
. As a result, in the general case, we obtain the relation referred below:

a = 1 -d 2sc 1 -d 1sc , (2.4.17)
where for the down-conversion d 2sc ∈ (0, 1) and d 1sc = 0 and for the up-conversion d 1sc ∈ (0, 1) and d 2sc = 0. Furthermore, considering that the Split-Pi converter is an electrical circuit with no interior power loss, we may apply a power conservation rule. Hence, the input power will be equal to the output power as referred below (see also Fig. 2.2.5):

P sc in =P sc out (2.4.18) v sc in i sc in =v sc out i sc out (2.4.19) i sc in i sc out = v sc out v sc in = α, (2.4.20)
where v sc in , by considering the resistance R 1sc added in the Bond graph (Fig. 2.4.7), is equal to v sc in = v si s R 1sc , respecting Ohm's law, and

v sc out = q 3sc C 3sc = v R Lsc . The v s , i s , q 3sc and v R Lsc can be found from the converter's PH model in (2.4.10)-(2.4.11).
For the simulation, we build the Split-Pi PH model in MATLAB/Simulink. We set the parameters according to Table 2.4.1 and we proceed to the simulations. In order to validate the PH model, we consider concrete values for the output voltage of the Split-Pi, u sc out , and we observe if, through the PH model, these values can be followed. Therefore, we consider, the first 300s, an output voltage equal to 13V , then we change the voltage to 15V and, finally, at 700s we change again to 14V in order to observe the response of the converter. Fig. 2.4.8 depicts the simulation results where the expected output voltage is well-followed by the Split-Pi PH representation. In Fig. 2.4.9, the corresponding duty cycle, d 2sc , is presented.

                          ṗ1cc ṗ2cc q1cc q2cc      = [J cc -R cc ]Q cc x cc +      1 0 0 0 0 0 0 1      v S -q2cc R Lcc C2cc , i S v R Lcc = 1 0 0 0 0 0 0 1 Q cc x cc , (2.4.21)
where

x cc = p 1cc p 2cc q 1cc q 2cc ∈ R 4×1 is the state vector, u cc = v S -q2cc R Lcc C2cc ∈ R 2×1
is the system's input vector denoting the input voltage v s from the DC bus and the load'

s current i R Lcc = q2cc R Lcc C2cc , y cc = i S v R Lcc ∈ R 2×1
is the output vector, i dc is the input current and v R Lcc is the output voltage of the converter. The matrix R cc is equal to 0 and J cc ∈ R 4×4 is depicted below:

J cc =     0 0 -(1 -d cc (t)) 0 0 0 -d cc (t) -1 1 -d cc (t) d cc (t) 0 0 0 1 0 0     .
(2.4.22)

The value d cc (t) ∈ (0, 1) is the control variable of the converter that corresponds to the duty cycle of the two converter's switches.

For the simulation, we construct the Ćuk converter PH model in MATLAB/Simulink. Based on Table 2.4.2, we move forward to the simulations. As for the Split-Pi converter, we consider several values for the expected output voltage v cc out and we observe if, through the PH model, they can be followed. Thus, we set an expected output voltage, firstly, equal to 8V , then changed to 12V and, finally, to 10V .

I 1cc , I 2cc 0.25 H C 1cc , C 2cc 0.0008 F R 0.025 Ω
For the simulation, similarly to (2.4.18) for the Split-Pi converter, we applied again the input/output power conservation for the Ćuk converter, as an electrical circuit without power losses, as follows:

P cc in =P cc out (2.4.23) v cc in i cc in =v cc out i cc out (2.4.24) i cc in i cc out = v cc out v cc in = d cc , ( 2

Central transmission network

This section analyzes the global central transmission network as a collection of meshgrids, which will be applied in the load balancing problem. This problem concerns the regulation of the DC-voltage in multiterminal DC networks and aims to achieve an optimal power flow (OPF) among the transmission lines, that is, finding the optimal working point of the system by using a set of constraints for the power, the current and the voltage [START_REF] Gavriluta | Hierarchical control of hv-mtdc systems with droop-based primary and opf-based secondary[END_REF]. In order to optimize the power flow of the system, the power losses in the transmission-line network, caused by the Joule effect, must be reduced , i.e., an associated cost function has to be minimized. Regarding the meshed DC microgrid considered in this work, the central transmission network is a cyclic network, because when the connecting nodes3 are linked together, they form cycles and, consequently, loops. In general, the transmission network is represented as a circuit with storage and dissipative elements. In here, the transmission lines, as already referred, are replaced by RL circuits as in Fig. 2.4.13. Every connecting node in this network can be connected to another meshgrid (Fig. 2.4.13) which leads to a meshed topology.

In order to analyze the network, we will develop a general PH state-space representation that can be used for each meshgrid separately. Then, the Bond graph of Fig. 2.4.1 will be used as an example to apply the

                                                 ṗ1c ṗ2c ṗ3c . . . ṗnl         = [J c -R c ]Q c x c + G c         v 1 v 2 v 3 . . . v nc         ,         i 1 i 2 i 3 . . . i nc         = G c Q c x c , (2.4.26)
where

x c = p 1c p 2c p 3c ... p nl ∈ R n×1 is the state vector, u c = v 1 v 2 v 3 ... v nc ∈ R node×1
is the input vector and y c = i 1 i 2 i 3 ... i nc ∈ R node×1 is the output vector of the system. The skewsymmetric matrix J c is equal to 0 for each meshgrid. As aforementioned, the J matrix (2.3.1) represents the oscillations among the electrical and the magnetic field. However, in this case, we consider only inductors, which lead to J c = 0. As regards the dissipation matrix, R c ∈ R n×4 , it is equal to:

R c = diag(R 1c , R 2c , R 3c , ..., R nl ), ( 2
                                                                     ṗ1 ṗ2 ṗ3 ṗ4 ṗ5 ṗ6 ṗ7             = [J c -R c ]Q c x c +             0 0 -1 1 0 0 0 0 -1 0 1 0 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 1 0 1 0 0 0 0 1 1 0                         v sn1 -v ug -v ev -v loads v pv -v b -v sn2             ,             i sn1 i ug i ev i loads i pv i b i sn2             =             0 0 1 1 0 0 0 0 -1 0 1 0 0 0 -1 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 -1 0 1 0 0 0 0 1 1 0             Q c x c , (2.4.29)
where

x c = p 1 p 2 p 3 p 4 p 5 p 6 p 7 ∈ R 7×1 is the state vector, u c = v sn1 -v ug -v ev - v loads v pv -v b -v sn2 ∈ R 7×1 is the input vector and y c = i sn1 i ug i ev i loads i pv i b i sn2 ∈ R 7×1
is the output vector of the system. In addition, the dissipation matrix R c ∈ R 7×7 is equal to:

R c = diag(R 1 , R 2 , R 3 , R 4 , R 5 , R 6 , R 7 ). (2.4.30)
Hence, the G c matrix contains all the power-preserving interconnections within the transmission line network and the R c matrix includes the power losses.

In the upcoming section, more details for the central transmission network will be given through a simplified architecture. The reduction of a model is usually considered as a forward step because of the systems' complexity. Sometimes the structure of a system has to be simplified [START_REF] Parisio | An mpc-based energy management system for multiple residential microgrids[END_REF][START_REF] Wang | Intelligent dc microgrid with smart grid communications: Control strategy consideration and design[END_REF] and replaced by a simpler dynamics, manageable from the viewpoint of further control design procedures and analysis.

Reduced model of the DC microgrid

In the sequel, the simplified version of the meshed DC microgrid will be described. To implement this simplification of the system, the number of the transmission lines will be reduced. This new scheme will facilitate the construction of the optimization problem in Chapter 4. The purpose is to consider a macroscopic view of the same meshed DC microgrid (Fig. 2.2.1), on which we will be able later to add the erased components. To proceed with this approach, the following assumptions have been made: reduction of the meshgrids down to one central transmission network; reduction of the transmission lines, setting R 5 = R 6 = R 7 equal to 0 and erasing I 5 = I 6 = I 7 . Moreover, considering the EVs station as the main energy storage, the batteries' section of the corresponding meshgrid can be removed and only the PV system remains; considering the EVs station as a storage system will account both for power storage and charging the EVs; Applying all these simplifications, we arrive to the system presented in Fig. 2.5.1. The meshed DC microgrid, at this point, contains a four-line transmission network with four nodes, on which the three power sources (UG, PV and ES) and the loads are connected. Below, we will deduce explicitly the dynamical models of the ES system (the battery linked together with the Split-Pi converter) and the latter proposed central transmission network.

Energy storage system

This subsection presents the mathematical model of the ES system which is highly relevant for the rest of the analysis and is required in order to provide a backup energy source. Below, the detailed electrical circuit of the KiBaM battery connected to the Split-Pi converter is introduced. In Fig. 2.5.2, we present the complete circuit of the ES system and its related Bond graph (Fig. 2.5.3). From the associated Bond graph, the PH formalism of the system is provided below:

ẋes (t) = [J es (d(t)) -R es ]Q es x es (t) + G es u es (t), y es (t) = G es Q es x es (t) + D es u es (t), (2.5.1)
where

x es (t) = p 1sc (t) p 2sc (t) q 1sc (t) q 2sc (t) q 3sc (t) q 1b (t) q 2b (t) ∈ R 7×1 , u es (t) = [-v DC (t) -i R 1b (t)] ∈ R 2×1 denotes the input voltage coming from the central network, otherwise called the DC bus of the DC microgrid, y es (t) = i DC (t) v R 1b (t) ∈ R 1×2
, where i DC (t) is the current during charging mode, and d(t) is the duty cycle of the converter's switches. Additionally, the diagonal matrix

Q es is equal to diag( 1/I 1sc , 1/I 1sc , 1/C 1sc , 1/C 2sc , 1/C 3sc , 1/C 1b , 1/C 2b ) ∈ R 7×7 .The skew-symmetric matrix J es (t) ∈ R 7×7 , the dissipation matrix R es ∈ R 7×7 , the G es ∈ R 2×7 and D es ∈ R 2×2 are equal to: J es (t) =           0 0 1 -(1 -d 1sc (t)) 0 0 0 0 0 0 (1 -d 2sc (t)) -1 0 0 -1 0 0 0 0 0 0 1 -d 1sc (t) -(1 -d 2sc (t)) 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0           , (2.5.2a)
R es =           0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 R1 sc 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 R 1b + 1 R 2b -1 R 2b 0 0 0 0 0 -1 R 2b 1 R 2b           , (2.5.2b) G es = 0 0 -1 R1 sc 0 0 0 0 0 0 0 0 1 -1 0 , D es = 1 R1sc 0 0 0 (2.5.2c)
where the duty cycles, d 1sc (t), d 2sc (t), are the control variables of the system. of charging, which is the summation of the available charge, q 1b , and the bound charge, q 2b , as in (2.5.1).

Note that for the simulation, we used the parameters illustrated in

Reduced transmission line network

Transmission lines modeled with RL circuits

Afterwards, we analyze the four-line reduced transmission network depicted in Fig. 2.5.1. To develop its PH model, we start from the PH model of the complete DC microgrid architecture (2.4.29) in order to produce the simplified version under the assumptions we made. Taking into account that there is no power loss in the transmission lines 5, 6, 7, the corresponding RL circuits will be erased. In addition, since only the part of the EV's will be regarded as an ES system, the v b variable will be equal to 0. The power resulting on the connecting node n = 5 and, subsequently, on the connecting node n = 4 will be generated by the PV system (v sn1 = v sn2 = v pv ). Therefore, considering that ṗ5sc , ṗ6sc , ṗ7sc and R 5 , R 6 , R 7 are equal to 0 in (2.4.29), we arrive at the following PH state-space representation:

                                  ṗ1 ṗ2 ṗ3 ṗ4      = [J c -R c ]Q c x c +      0 0 -1 1 0 -1 0 1 1 0 1 0 1 1 0 0           v pv -v ug -v es -v loads      ,      i pv i ug i es i R L      =      0 0 1 1 0 -1 0 1 -1 0 1 0 1 1 0 0      Q c x c , (2.5.3)
where

x c = p 1 p 2 p 3 p 4 ∈ R 4×1 is the state vector, u c = v pv -v ug -v es -v loads ∈ R 4×1 is
the input vector and y c = i pv i ug i es i loads ∈ R 4×1 is the output vector of the system. The skewsymmetric matrix J c is equal to 0 and the dissipation matrix R c ∈ R 4×4 is equal to: where for linear maps F :F → V and E:E → V obeying to:

R c = diag(R 1 , R 2 , R 3 , R 4 ). ( 2
EF * + F E * = 0, rank(F + E) = dimF, (2.5.6)
where V is a linear space of equal dimension with F, and where F * :V * → E and E * :V * → (F * ) * = F represent the adjoint maps of F and E respectively. The linear maps F and E are described by n L × n L matrices F and E which satisfy:

EF + F E = 0, rank(F + E) = dimF,
which is the matrix kernel representation.

For the cyclic reduced network (Fig. 2.5.6), the kernel representation is given below:

                    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 -1 0 0 0 0 0 0 0 0 1 0 0 0 -1 0 0 0 0 0 0 0 0 1 0 0 0 -1 0 0 0 0 0 0 0 0 1 0 0 0 -1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 -1 0 -1 0 1 0 0 0 0 0 1 0 0 0 0 -1 -1 0 0 0 0 1 0 0 0 -1 0 1 0 0 0 0 0 0 0 -1 0                     F                     i 1 i 2 i 3 i 4 i R1 i R2 i R3 i R4 i pv i ug i es i loads                     f + +                     -1 0 0 0 -1 0 0 0 0 0 1 -1 0 -1 0 0 0 -1 0 0 0 -1 0 -1 0 0 -1 0 0 0 -1 0 -1 0 -1 0 0 0 0 -1 0 0 0 -1 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                     E                     v 1 v 2 v 3 v 4 v R1 v R2 v R3 v R4 -v pv -v ug v es v loads                     e =                     0 0 0 0 0 0 0 0 0 0 0 0                     . ( 2
D = {(f, e) ∈ F × E|f = Je + Gλ, G = 0, λ ∈ V} , (2.5.8)
with a skew-symmetric mapping J : F → E, a linear mapping G such that imG = {f ∈ F|(f, 0) ∈ D} and kerJ = {e ∈ E|(0, e) ∈ D}. Assuming that F = m L (≤ n L ), then m L independent columns of F can be selected and placed in a matrix

F 1 . As a next step, F = [F 1 | F 2 ] and, accordingly, E = [E 1 | E 2 ]
, where f = f 1 f 2 and e = e 1 e 2 . Hence, the matrix [F 1 | E 2 ] is invertible and the Dirac structure can be written as follows:

D = f 1 f 2 ∈ F, e 1 e 2 ∈ E| f 1 e 2 = J e 1 f 2 , ( 2
.5.9)

where

J:=-[F 1 | E 2 ] -1 [F 2 | E 1 ]
which is a skew-symmetric matrix.

                    i R1 i R2 i R3 i R4 i pv i ug i es i loads v 1 v 2 v 3 v 4                     =                     0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 -1 -1 0 0 0 0 0 0 0 0 0 -1 0 1 0 0 0 0 0 0 0 0 1 0 -1 0 0 0 0 0 0 0 0 0 -1 -1 0 0 -1 0 0 0 0 0 -1 1 0 0 0 0 0 -1 0 0 0 1 0 1 0 0 0 0 0 0 -1 0 1 0 1 0 0 0 0 0 0 0 0 -1 1 -1 0 0 0 0 0 0                                         v R1 v R2 v R3 v R4 -v pv -v ug v es v loads i 1 i 2 i 3 i 4                     .
(2.5.10)

Additionally to the PH model, we always take into consideration the power conservation equation of the transmission network written below:

P ug + P pv -P es -P loads -P R1 -P I1 -P R2 -P I2 -P R3 -P I3 -P R4 -P I4 = 0. (2.5.11)
The problem in this form becomes complicated and many difficulties can still emerge during the construction of the optimization problem. Therefore, eliminating a number of states in the transmission network would be useful. In the subsequent section, a further reduction of the transmission-line network will be proposed.

Replacing transmission lines by resistors

This section presents a further simplification of the central transmission network according to which the inductors, I 1 , I 2 , I 3 and I 4 will be eliminated. In DC networks, the existence of the inductors is necessary only when the current flow of the system starts increasing. But once the system comes to its steady state, where the current flow is stable, the inductors have no effect in the DC networks. Furthermore, in DC networks the distances among the components are small meaning that the length of the transmission lines is negligible. Together with the continuous development of the DC equipment, these are two factors which lead towards the decrease of the power loss within DC networks [START_REF] Elsayed | Dc microgrids and distribution systems: An overview[END_REF][START_REF] Lotfi | Ac versus dc microgrid planning[END_REF].

Therefore, in this work, to study the load balancing problem of the central transmission network of the simplified DC microgrid (Fig. 2.5.1), the inductors will be discarded. Thus, we illustrate the correspondinglymodified Bond graph for the central DC bus in Fig. 2.5.7. Hence, from the kernel representation in (2.5.7), erasing the inductors' variables (current, i 1 , i 2 , i 3 , i 4 , and voltage, v 1 , v 2 , v 3 , v 4 ), we obtain the following for the Bond graph in Fig. 2.5.7:

            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 -1 0 -1 0 1 0 1 0 0 0 0 -1 -1 1 0 0 0 -1 0 1 0 0 0 -1 0             F             i R1 i R2 i R3 i R4 i pv i ug i es i loads             f +             -1 0 0 0 0 0 1 -1 0 -1 0 0 0 -1 0 -1 0 0 -1 0 -1 0 -1 0 0 0 0 -1 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0             E             v R1 v R2 v R3 v R4 -v pv -v ug v es v loads             e =             0 0 0 0 0 0 0 0             .
(2.5.12) Similarly, from (2.5.12), the hybrid input-output representation of the cyclic network in Fig. 2.5.7 is:

            i pv i ug i es i loads v R1 v R2 v R3 v R4             =             0 0 0 0 0 0 1 1 0 0 0 0 0 1 0 -1 0 0 0 0 -1 0 1 0 0 0 0 0 1 1 0 0 0 0 1 -1 0 0 0 0 0 -1 0 -1 0 0 0 0 -1 0 -1 0 0 0 0 0 -1 1 0 0 0 0 0 0                         -v pv -v ug v es v loads i R1 i R2 i R3 i R4             .
(2.5.13) Furthermore, the power conservation equation of the transmission network (Fig. 2.5.7) becomes:

P ug + P pv -P es -P loads -P R1 -P R2 -P R3 -P R4 = 0.

(2.5.14)

From the Bond graph in Fig. 2.5.7, the PH representation is presented below, where there are no states and only dissipative elements are included:

    i pv i ug i es i loads     =           1 R 3 + 1 R 4 - 1 R 4 - 1 R 3 0 1 R 4 - 1 R 4 - 1 R 2 0 1 R 2 1 R 3 0 - 1 R 3 - 1 R 1 1 R 1 0 1 R 2 1 R 1 - 1 R 1 - 1 R 2               v pv v ug v es v loads     .
(2.5.15) Therefore, the following equations describe for the current flows: from node 1:

i loads = i R1 + i R2 = v es -v loads R 1 + v ug -v loads R 2 ;
(2.5.16) from node 2:

i es = i R3 -i R1 = v pv -v es R 3 - v es -v loads R 1 ;
(2.5.17) from node 3:

i ug = i R4 -i R2 = v pv -v ug R 4 - v ug -v loads R 2 ;
(2.5.18) from node 4:

i pv = i R3 -i R4 = v pv (t) -v es (t) R 3 - v ug (t) -v pv (t) R 4 . (2.5.19)
Through (2.5.13), (2.5.16), (2.5.17), (2.5.18), (2.5.19), every power variable of the DC network can be expressed below in function of the voltages on the connecting nodes:

P ug = v ug i ug = v ug [i R4 -i R2 ] = v ug v ug -v pv R 4 - v loads -v ug R 2
(2.5.20a)

P es = v es i es = v es [i R3 -i R1 ] = v es v pv -v es R 3 - v es -v loads R 1
(2.5.20b)

P pv = v pv i pv = v pv [i R3 -i R4 ] = v pv v pv -v es R 3 - v ug -v pv R 4
(2.5.20c)

P loads = v loads i loads = v loads [i R1 -i R2 ] = v loads v es -v loads R 1 - v loads -v ug R 2 , (2.5.20d)
where P ug , P es , P pv , P loads are the powers of the UG, ES, PV and loads systems respectively. Furthermore, regarding the power losses within the DC bus, P R1 , P R2 , P R3 , P R4 , the following relations are also deduced:

P R1 = [v es -v loads ] 2 R 1 , (2.5.16e
)

P R3 = [v pv -v es ] 2 R 3 , (2.5.16f
)

P R2 = [v loads -v ug ] 2 R 2 , ( 2
.5.16g)

P R4 = [v ug -v pv ] 2 R 4 . (2.5.16h)
The previous relations (2.5.20a-2.5.16h) are important and will be applied in the construction of the objective function for the power loss minimization problem in chapter 4.

Optimization objectives and constraints

Before starting to elaborate the optimization problem, it is important to present the overall electrical circuit of the simplified version of the DC microgrid (Fig. 2.5.1). Fig. 2.6.1 introduces all the inner interconnections of the circuits of the sources and the Split-Pi converters.

Objectives

The main objective of this work is to solve the load balancing problem. To achieve this, the problem will be divided into two sub-problems: the cost minimization, according to which the electricity cost of the UG power purchase will be penalized. The goal is to sell power to the UG, generated by the renewable resources, and to exploit the ES system towards the consumers' benefit. The cost function which penalizes the electricity cost is written below:

J 1 (d(t)) = t f t0 ep(t)P ug = t f t0 ep(t)(P loads + P es -P pv ) = t f t0 ep(t)(P loads + v es i es -P pv ), (2.6.1)
where reference profiles will be taken into account for the PV and the loads. The P ug is replaced by the power conservation among the sources and loads (P ug = P loads + P es -P pv ), without considering the power loss (2.5.14); the optimal power flow, minimizing the power dissipation in the DC bus. This problem aims to find the best possible direction for the power to flow in order to have the less possible power losses during transmission. Therefore, the cost function which minimizes the energy dissipation will be the following:

J 2 (d(t)) = - t f t0 P R = - t f t0 (P R1 + P R2 + P R3 + P R4 ), (2.6.2)
where P R1 , P R2 , P R3 , P R4 corresponds to the power loss inside the central network (2.5.16e-2.5.16h), for which again the PH model through hybrid input-output representation of the central transmission network will be considered. Relation (2.5.14) is taken into account.

We observe that in both cases the power conservation is taken into account and the control variable is the duty cycle of the switches in the converters, d(t). The ES system and the central transmission network will be analytically included in the optimization problem, utilizing the PH models previously developed in this chapter.

Constraints

This subsection will describe the general constraints which will be taken into account for the energy management problem. As already mentioned, the ES system is a major component of the overall scheme and will be thoroughly analyzed. The batteries have a limited lifetime which demands to respect certain constraints concerning the battery's characteristics. Therefore, the constraints considered in this work for the ES system are (see also Fig. 2.6.1):

q min 1b ≤ q 1b (t) ≤ q max 1b , q min 2b ≤ q 2b (t) ≤ q max 2b , Charging mode : v min b,charging ≤ v b,charging (t) ≤ v max b,charging , i min b,charging ≤ i b,charging (t) ≤ i max b,charging , Dischaging mode : v min b,dicharging ≤ v b,dicharging (t) ≤ v max b,dicharging , i min b,dicharging ≤ i b,dicharging (t) ≤ i max b,dicharging .
In addition to these constraints, the output voltage on the connecting nodes in the DC-bus must be always close to a certain voltage value, 400 V in this case, because we consider DC networks. Therefore, according to Fig. 2.6.2, the limitations considered for the connecting nodes are given below:

v min,h DC ≤v ug (t) ≤ v max,h DC ,
(2.6.3a)

v min,h DC ≤v pv (t) ≤ v max,h DC , (2.6.3b) v min,h DC ≤v es (t) ≤ v max,h DC , (2.6.3c) v min,h DC ≤v loads (t) ≤ v max,h DC .
(2.6.3d) (2.6.5)

Finally, the external grid cannot consume an internal amount of power. An upper and a lower limit must be defined as follows:

P min ug ≤ P ug (t) ≤ P max ug , (2.6.6)

with P ug (t) = v ug (t)i ug (t).

Conclusions

In this chapter, the meshed DC microgrid architecture of the project was presented. Afterwards, we introduced the modeling methodology used to represent the components of the system. The Bond graph and port-Hamiltonian representations were developed proceeding to the detailed analysis of the system as a collection of electrical circuits. The global dynamics of the system was also formulated by focusing on the central transmission network.

To simplify the complexity of the system for the upcoming analysis, we established a simplified architecture of the proposed meshed DC microgrid reducing the number of transmission lines. Studying each meshgrid separately leads to a simpler PH model for the power transmission among the components. Through this representation, the explicit dynamical model of the central transmission network and the ES system are presented, important for the development of the multi-scale optimization problem.

The highlights of this chapter are enumerated below: development of a meshed DC microgrid using Split-Pi converters for voltage regulation and power flow direction in the network; development of the dynamical models using the Bond graphs and their associated port-Hamiltonian representations for each component of the DC microgrid. The dynamics for each component of the system has been introduced explicitly. A similar approach was followed for the central transmission network, thus providing the analytical framework of the power line-transmission;

suitable cost functions and constraints have also been presented.

The next chapter will investigate the inverse dynamics of the PH models developed here. This is an important step which will provide directly the states of the system, maintaining its properties, to be used in the optimization problem. Towards this direction, two methods will be considered, differential flatness and bicausality.

Chapter 3

From port-Hamiltonian to differential flatness representation

Introduction

This chapter presents the flat representation of the PV, ES and loads components of the reduced DC microgrid model (Fig. 2.5.1) detailed in chapter 2. Differential flatness was first introduced by Michel Fliess [START_REF] Fliess | Differential flatness and defect: an overview[END_REF] and his co-researchers in the late 1980s. It is an advantageous method, which can efficiently deal with control and optimization problems of nonlinear dynamical systems [START_REF] Levine | Analysis and control of nonlinear systems: A flatness-based approach[END_REF][START_REF] Murray | Differential flatness of mechanical control systems: A catalog of prototype systems[END_REF]Rigatos, 2015a).

There are several works in the literature in different domains that employ the differential flatness tracking control design. Many researchers use its properties in motion planning problems (Hervagault et As already mentioned in chapter 2, through port-Hamiltonian (PH) formulation we obtain an explicitly described and well-structured state-space representation of the system with a set of differential equations. To proceed to its supervision, it is necessary to be able to extract several physical quantities such as the voltage, the current, the power and the like. Therefore, differential flatness can be considered as a suitable tool to inverse the system dynamics generating in this way all the states and inputs in function of the flat outputs of the system and a finite number of their derivatives. It is a method closely associated to the notion of controllability, since according to Fliess et al. (1995), a nonlinear system is flat if and only if it is controllable. Due to its properties, it can be easily combined with several methods of parametrization such as the B-spline parametrization, which accounts for constraint validation in continuous-time.

Another important issue, that concerns the researchers, is the calculation of the flat outputs which can be very complex and not straightforward to find. In the literature, many useful approaches are provided, such as the algorithm for flat representation proposed in [START_REF] Franke | On the computation of flat outputs for nonlinear control systems[END_REF]. The algorithm provides an analytical computation of flat outputs for nonlinear control systems and will be described and implemented later in this chapter. Apart from the algorithm, another method was introduced in Gil et al. (1997) and [START_REF] Richard | Analysis of flatness using bond graphs and bicausality[END_REF], which proves the flatness of a system after inversing the dynamics of its associated Bond graph. This method is called bicausality, firstly presented by [START_REF] Gawthrop | Bicausal bond graphs[END_REF]. Bicausality is a supplementary approach to the existing theory of the Bond graphs which inverses the system dynamics by decoupling the pairs of efforts and flows at each bond. Hence, bicausality could be an appropriate method to investigate flatness and, subsequently, develop a more straightforward way to find the possible flat outputs of the PH systems.

In this chapter, we will proceed by briefly introducing the differential flatness theory. Afterwards, we will apply the theory on the reduced DC microgrid system (section 2.5) using two sets of flat outputs, one derived from the algorithm proposed in Franke and Robenack (2013) and another considering a random flat output set. The proposed sets of flat outputs will be validated through simulations in MATLAB/Simulink. Next, we will present the bicausal Bond graph and we will demonstrate a method which can result in the generation of the possible flat outputs for a particular class of PH systems. The method will be explained through illustrative examples of electrical circuits. Fig. 3.1.1 presents the methods used in this chapter which both results in the inverse dynamics of the system presenting the states and inputs in function of the flat outputs and their derivatives.

Differential flatness

Differential flatness is a structural property of a class of nonlinear dynamical systems, denoting that all system variables (the states and the control inputs) can be written in terms of a set of specific variables, the so-called flat outputs (equal in number to the number of inputs), and their derivatives (Rigatos (2015b)).

ẋ(t) = f (x(t), u(t)), (3.2.1)
where the state vector x(t) ∈ R n and the input vector u(t) ∈ R m . In addition, f (0, 0) = 0 and rank ∂f ∂u = m must be verified. The system can be characterized as differentially flat, if there exists a flat output vector:

z(t) = [z 1 (t) z 2 (t) ... z m (t)] (3.2.2)
which satisfies the conditions below:

the flat output z(t) is represented in function of the states and the inputs of the system and their derivatives, z(t) = Φ(x(t), u(t), u(t), ü(t), ..., u (k+1) (t));

the states and the inputs of the system are described in terms of the flat outputs and a finite number of their derivatives, x(t) = Φ x (z(t), ż(t), z(t), ..., z (k) (t)) and u(t) = Φ u (z(t), ż(t), z(t), ..., z (k+1) (t));

the flat outputs z(t) and their derivatives are differentially independent, which means that they cannot satisfy equations of the form Φ(z(t), ż(t), z(t), ..., z (k+1) (t)) = 0.

Differential flatness has many advantages due to its properties, two of them, useful for this work, are presented below: the number of the control inputs defines the number of the flat outputs of the system ((Lévine, 2011)); it is commonly established that a flat system is considered as linearizable through a dynamic feedback and an appropriate change of coordinates. In general, when a system is linear and controllable, it is, consequently, flat. Therefore, a nonlinear system that can be converted into a linear, controllable system is also flat. Note that the previous statements do not prove that flatness is a method of linearization. However, when a nonlinear system is flat, its structure can be inscribed in a way to be employed in control for trajectory generation, planning, stabilization etc. [START_REF] Levine | Analysis and control of nonlinear systems: A flatness-based approach[END_REF][START_REF] Rigatos | Nonlinear control and filtering using differential flatness approaches: applications to electromechanical systems[END_REF] Proof. Consider the general state-space representation of a PH system:

ẋ(t) = [J(t) -R]Qx(t) + Gu(t), (3.2
.3a)

y(t) = G Qx(t) + Du(t), (3.2.3b)
where x(t) ∈ R n is the state vector, u(t) ∈ R m is the input vector and y(t) ∈ R m is the output vector. Moreover, J(t) ∈ R n×n is the interconnection matrix, R ∈ R n×n is the dissipative matrix, G ∈ R n×m is the control matrix and D ∈ R m×m is the input matrix. For the Hamiltonian, H, and the J(t) -R matrix the following are considered:

the Hamiltonian,

H(x) = 1 2 x(t) Qx(t) (2.3.2)
, is in quadratic form and positive definite, since Q is positive definite and invertible. Its gradient vector ∇H = ∂ x H(x) will be ∂ x H(x) = Qx(t); J(t) -R is square and invertible.

Considering as flat outputs a number of the states of the system, equal in number to the number of control inputs (see also Definition 3.2.1), the flat outputs can be written as z M (x) = Φ(x 1 (t), x 2 (t), ..., x n (t)), where M = 1, 2, .., m, and the states can be written as x N (z) = Φ(z 1 (t), z 2 (t), ..., z n (t), ż1 (t), ż2 (t), ..., żn (t), ..., z The matrix Q M ∈ R n×n , n = m is also square and invertible. Then, the control inputs and the states can be written in function of the flat outputs and their derivatives as follows:

u M (z) = G -1 M ( ẋN (z) -[J(t) -R]Q M x N (z)), (3.2.4) x N (z) = [J(t) -R] -1 Q -1 ( ẋN (z) -Gu M (z)). (3.2.5)
Additionally, referring to the controllability of the PH systems considered here, the following remark is added:

Remark 3.2.3. (Port-controlled Hamiltonian system (Maschke and van der Schaft, 1992)): The interaction of a PH system with its environment is described defining external ports in the generalized Dirac structure (DS), similarly to the definition of the energy storing elements (see also Theorem 2.3.3). The control inputs can be represented as generalized effort or flow sources connected to the external ports of the model, with its effort or flow variables being the inputs. An external control has to act through the ports of the system, in order to change its states and, hence, its energy. Therefore, the system is port-controlled. In the network's framework an output is naturally generated, associated with each port-input control effort or flow, namely the conjugated flow or effort on the port respectively.

Consequently, from Definition 3.2.1, Proposition 3.2.2 and Remark 3.2.3, the PH systems, considered in this work, are controllable and flat.

Flat representation of the reduced DC microgrid system

Next, the differential flatness of the different components (ES, PV, loads) of the reduced DC microgrid system (section 2.5.1) will be studied taking into account their state-space representations from chapter 2. Two sets of flat outputs for each part will be considered, one set generated by the algorithm proposed in [START_REF] Franke | On the computation of flat outputs for nonlinear control systems[END_REF] and one set by randomly choosing a number of the system's states according to Proposition 3.2.2. In the following, the steps of the algorithm will be described and, afterwards, the flat representations of the ES, the PV and the loads systems will be analyzed. We highlight that the flat representation of the system will be used afterwards in chapter 4 for the optimal profile generation.

Algorithm description for flat output representation

The algorithm, implementing a symbolic linearization of the system and following a consecutive calculation of matrices, nullspaces and inverses, concludes to the corresponding sets of flat outputs of the system through integration. To begin, it considers control systems which are described implicitly through the implicit function theorem (Krantz and Parks, 2012) referred below. 2 , ..., x 0 n , y 0 1 , y 0 2 , ..., y 0 m )∈ R n × R m and on its neighborhood differentiable functions F 1 , F 2 , ..., F n . If F 1 (x 0 , y 0 ) = F 2 (x 0 , y 0 ) = . . . = F n (x 0 , y 0 ) = 0 and the n × n matrix at (x 0 , y 0 ):

          ∂F 1 ∂x 1 ∂F 1 ∂x 2 . . . ∂F 1 ∂x n ∂F 2 ∂x 1 ∂F 2 ∂x 2 . . . ∂F 2 ∂x n . . . . . . . . . . . . ∂F n ∂x 1 ∂F n ∂x 2 . . . ∂F n ∂x n           (3.3.1)
is different from 0, then there is a neighborhood U of y 0 = (y 0 1 , y 0 2 , . . . , y 0 m ) ∈ R m , a neighborhood V of x 0 = (x 0 1 , x 0 2 , . . . , x 0 n ) ∈ R n and a unique mapping ϕ : U -→ V for which ϕ(y 0 ) = x 0 and F 1 (ϕ(y), y) = F 2 (ϕ(y), y) = ... = F n (ϕ(y), y) = 0 for all y ∈ U , where ϕ is differentiable. Definition 3.3.2. (Implicit control systems [START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF]): For a n-manifold X, which is infinite and differentiable, its tangent space at a point x ∈ X is expressed by T x X and its tangent bundle T X. For a function F , which belongs to a set of C ∞ mappings (C ∞ ∈ (T X; R n-m ), the implicit system is considered as:

F (x, ẋ) = 0, (3.3.2)
for which the rank( ∂F ∂ ẋ ) = nm in a suitable open subset of T X. Remark 3.3.3. According to Theorem 3.3.1, any explicit system as in (3.2.1), where f (x, u) ∈ T x X for every x ∈ X and u in an open subset U of R m , and rank( ∂f ∂u ) = m in a suitable subspace of X × U , can be locally converted to (3.3.2) and vice versa [START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF]).

In the following, the steps of the algorithm for flat representation, as used in this work, will be presented according to Franke and Robenack (2013):

Algorithm for flat representation 1. The implicit system is introduced as in (3.3.2), based on Definition 3.3.2 and Remark 3.3.3:

F (x, ẋ) = 0, x ∈ R n , (3.3.3)
of nm equations, where n is the number of states and m is the number of inputs. The implicit control system is obtained after the elimination of the inputs of the system (Lévine, 2004);

2. Then, a dynamic feedback linearization is applied to the system (3.3.3) as follows:

n i=1 ∂F ∂x i dx i + ∂F ∂ ẋi d ẋi = dF = 0, (3.3.4)
from which the corresponding tangent system is found. According to (Aranda-Bricaire et al., 1995; Franke and Robenack, 2013), the tangent system consists of a number of one-forms ω equal in number to the number of the inputs of the system, m:

ω j = n i=1 F i (x)dx i , j = 1, 2, . . . , m. (3.3.5)
Each one-form, ω j , is called tangent flat output of the tangent system. The goal of the proposed algorithm is to find a possible set of tangent flat outputs from which the coordinates dx can be computed through integration. In general, the algorithm for flat representation is based on an iteration process that eliminates gradually the dimensions of the tangent system proving at the same time its controllability and achieving eventually to generate the set of the tangent flat outputs.

3. For the sake of convenience, equation (3.3.4) is replaced by (nm) × n matrices, P 0,[i] , P 1,[i] , and n × 1 matrices, u [i] , u[i] , where:

u [i]
contains the dx i coordinates and u[i] the d ẋi coordinates of the tangent system;

P 0,[i] is composed by ∂F ∂x i (x i , ẋi ) and P 1,[i] is composed by ∂F ∂ ẋi (x i , ẋi ) derivatives.
Consequently, (3.3.4) is rewritten as follows:

P 0,[i] u [i] + P 1,[i] u[i] = 0, (3.3.6)
where

u [i] = P + 1,[i] u [i+1] + P ⊥ 1,[i] w [i+1] , (3.3.7
)

u[i] = P + 1,[i] u[i+1] + Ṗ + 1,[i] u [i+1] + P ⊥ 1,[i] ẇ[i+1] + Ṗ ⊥ 1,[i] w [i+1] . (3.3.8)
Note that i is the step of the iteration process within the algorithm. P + 1,[i] matrix is the right inverse of P 1,[i] and P ⊥ 1,[i] matrix is the orthogonal complement of P 1,[i] which confirm: , it obeys to the statements written below:

P 1,[i] P + 1,[i] = I, P 1,[i] P ⊥ 1,[i] = 0. ( 3
u[i+1] + P 0,[i] -Ṗ1,[i] P + 1,[i] A [i] u [i+1] + P 0,[i] -Ṗ1,[i] P ⊥ 1,[i] B [i] w [i+1] = 0. ( 3
u [i] = P + 1,[i] u [i+1] + P ⊥ 1,[i] w [i+1] + Z [i] z [i+1] . ( 3
(P 0,[i] -Ṗ1,[i] )Z [i] = 0 and P 1,[i] Z [i] = 0; the inverse of the Z [i] matrix, Z + [i] is equal to z [i+1] = Z + [i] u [i]
and verifies the following conditions:

Z + [i] Z [0] = I, Z + [i] P ⊥ 1,[i] = 0, Z + [i] P + 1,[i] = 0. ( 3
B ⊥ [i] A [i] u [i+1] + B ⊥ [i] u[i+1] = 0, (3.3.13)
which results in the next step i + 1 of the iteration process where the dimensions of (3.3.6) are reduced. As a consequence, (3.3.6) is replaced by: :

P 0,[i+1] u [i+1] + P 1,[i+1] u[i+1] = 0. ( 3
w [i+1] = -B + [i] A [i] u [i+1] -B + [i] u[i+1] . ( 3
u [i+1] = P 1,[i] u [i] . (3.3.16)
The previously mentioned process continues until the matrix B [i] has full column rank, which concludes to a possible set of flat outputs, ω j of the tangent system:

ω j = P 1,[i] Z + [i] u [i] . (3.3.17)
Consequently, integrating (3.3.17) leads to a possible set of flat outputs for the system ẋ(t) = f (x(t), u(t)) (3.2.1).

Flat representation of the ES system

In this section, the ES model (the interconnected system of the battery together with the Split-Pi converter) will be considered (2.5.1, 2.5.2a, 2.5.2b, 2.5.2c). Its state-space representation is written below (Fig. 3.3.1), already presented in PH form in section 2.5.1: To proceed to the flat output calculation, we, primarily, observe that the ES system of the microgrid has:

ṗ1sc (t) = q 1sc (t) C 1sc - q 2sc (t) C 2sc (1 -d 1sc (t)), (3.3.18a) ṗ2sc (t) = q 2sc (t) C 2sc (1 -d 2sc (t)) - q 3sc (t) C 3sc , (3.3.18b)
q1sc (t) = v s R 1sc - q 1sc (t) C 2sc R 1sc - p 1sc (t) I 1sc , (3.3.18c)
q2sc (t) = p 1sc (t) I 1sc (1 -d 1sc (t)) - p 2sc (t) I 2sc (1 -d 2sc (t)), (3.3
.18d) q3sc (t) = p 2sc (t) I 2sc -i R 1b (t), (3.3.18e) q1b (t) = i R 1b (t) - q 1b (t) C 1b R 2b + q 2b (t) C 2b R 2b , (3.3.18f) q2b (t) = q 1b (t) C 1b R 2b - q 2b (t) C 2b R 2b . ( 3
seven states (3.3.18a-3.3.18g); p 1sc (t), p 2sc (t), q 1sc (t), q 2sc (t), q 3sc (t), q 1b (t), q 2b (t), four inputs (3.3.18a-3.3.18g):

d 1sc (t), d 2sc (t), v s (t), i R 1b (t),
and the circuit's parameters and the values of the resistors are considered as constants,

1/C 1sc = α, 1/C 2sc = β, 1/C 3sc = γ, 1/C 1b = δ, 1/C 2b = , 1/I 1sc = η, 1/I 2sc = µ, 1/R 1sc = ν, 1/R 1b = σ, 1/R 2b = ψ).
By taking into account the previously mentioned conditions, we derive the ensuing equations from the state-space representation (3.

3.18a)-(3.3.18g) 1 : ṗ1sc = αq 1sc -βq 2sc (1 -d 1sc ), (3.3.19a) ṗ2sc = βq 2sc (1 -d 2sc ) -γq 3sc , (3.3.19b) q1sc = νv s -ανq 1sc -ηp 1sc , (3.3.19c) q2sc = ηp 1sc (1 -d 1sc ) -µp 2sc (1 -d 2sc ), (3.3.19d) q3sc = µp 2sc -i R 1b , (3.3.19e) q1b = i R 1b -δψq 1b + ψq 2b , (3.3.19f) q2b = δψq 1b -ψq 2b . (3.3.19g)
Afterwards, we will elaborate two sets of flat outputs for the ES system, one generated by the algorithm for flat representation as presented in section 3.3.1 and one chosen randomly according to the Proof 3.2. Finally, we will compare the results through simulations.

βq 2sc q2sc + ηp 1sc ṗ1sc -αηp 1sc q 1sc + µp 2sc ṗ2sc + γµq 3sc p 2sc = 0, (3.3
.20a) q1b + δψq 1b + ψq 2b -µp 2sc + q3sc = 0, (3.3.20b)
q2bδψq 1b + ψ q 2b = 0. The detailed explanation of the algorithm for flat representation procedure for the ES system is presented in Appendix A.1.1. According to (A.1.2b) and (A.1.12), we conclude to the subsequent tangent flat outputs:

    ω 1 ω 2 ω 3 ω 4     = P 1,[0] Z + [0]           dp 1sc dp 2sc dq 1sc dq 2sc dq 3sc dq 1b dq 2b           , (3.3.22)
where From the algorithm, we conclude to a possible set of flat outputs through the tangent system of the ES. However, according to [START_REF] Levine | Analysis and control of nonlinear systems: A flatness-based approach[END_REF], multiple flat output sets may exist for each physical system. In the following section, we consider another set of flat outputs for the ES and we compare with the preceding results.

P 1,[0] =   ηp 1sc µp 2sc 0 βq 2sc 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1   , Z + [0] = 0 0 0 1 0 0 0 . q 1b = Φ 6 (z 3 , ż3 ), (3.3.27f)
q 2b = Φ 7 (z 3 ), (3.3.27g) d 1sc (t) = Φ 8 (z 1 , z 2 , ż2 , z2 , z 3 , ż3 , z3 , z 4 ), (3.3.27h) d 2sc (t) = Φ 9 (z 2 , z2 , z 3 , ż3 , z3 , z 4 ), (3.3.27i) v s = Φ 10 (z 1 , ż1 , z1 , z 2 , ż2 , z2 , ... z 2 , z 3 , ż3 , z3 , ... z 3 , z 4 , ż4 , z4 ), (3.3.27j) i R 1b = Φ 11 ( ż3 , z3 ). ( 3

ES flat representation through another set of flat outputs

In this section, taking into account the Proof 3.2, we suggest another set of flat outputs. According to [START_REF] Levine | Analysis and control of nonlinear systems: A flatness-based approach[END_REF], the number of control inputs equals the number of flat outputs which can be in function of the states of the system (Definition 3.2.1). Therefore, for the ES PH system of (3.3.18a)-(3.3.18g), we consider four states as flat outputs in function of which the system can be written. The selected set is the following:

z = p 1sc q 3sc q 2b q 2sc . (3.3.28)
Afterwards, considering the aforementioned set of flat outputs (3.3.28), the states and control inputs of the system are rewritten in flat representation as shown below:

p 1sc = z 1 , (3.3.29a
)

p 2sc = 1 µ ( ż2 + 1 δψ z3 + ( δ + 1) ż3 ), (3.3.29b
)

q 1sc = 1 a ż1 + β α z 4 ( ż4 + µ( 1 µ ( ż2 + 1 δψ z3 + ( δ + 1) ż3 )( γz 2 βz 4 + 1 µβ z2 z 4 + 1 µδψβ ... z 3 z 4 + 1 µβ ( δ + 1) z3 z 4
))

ηz 1 ),

(3.3.29c)

q 2sc = z 4 , (3.3.29d)
q 3sc = z 2 , (3.3.29e
)

q 1b = 1 δψ ż3 + δ z 3 , (3.3.29f
)

q 2b = z 3 , (3.3.29g
)

d 1sc = 1 - ż4 + µ( 1 µ ( ż2 + 1 δψ z3 + ( δ + 1) ż3 )( γz 2 bz 4 + 1 µβ z2 z 4 + 1 µδψβ ... z 3 z 4 + 1 µβ ( δ + 1) z3 z 4
))

ηz 1 , (3.3.29h
)

d 2sc = 1 - γz 2 βz 4 - 1 µβ z2 z 4 - 1 µδψβ ... z 3 z 4 - 1 µβ ( δ + 1) z3 z 4 , (3.3.29i) v s = 1 n q1sc + αq 1sc + η ν p 1sc , (3.3.29j) i R 1b = 1 δψ z3 + ( δ + 1) ż3 . (3.3.29k)
Moreover, the compact flat representation of the states and inputs is presented in function of the flat outputs (3.3.28) and their derivatives: Next, we proceed similarly to the simulations in order to validate the second set of flat outputs (3.3.28), considering the same reference values for the control inputs of the ES PH system, d 1sc , d 2sc , v s , i R 1b as before. Fig. 3.3.4 depicts the signal response of the control inputs which approaches very well the reference profiles. Consequently, both sets of flat outputs (3.3.25a-3.3.25d and 3.3.28) can be used to generate an equivalent flat system for the ES. However, the flat representation generated by the algorithm (3.3.26a-3.3.26l) seems to be a better choice, since in the control inputs simulations (Fig. 3

p 1sc = Φ 12 (z 1 ), ( 3

Flat outputs representation of the PV and the loads

In this section, the flat outputs for the PV and the loads' system (see also Fig. 2.5.1) will be investigated, in the same way as previously for the ES system, by computing two different flat output sets, one with the algorithm and another through a selection among the states of the system. The PV and the loads are two different systems composed by (see also Fig. 2.6.1): in the case of the PV, the PV is the power source which generates the input voltage, v s , of the Split-Pi converter (Fig. 3.3.5) and at the output there is the central transmission network, the common DC-bus. In Fig. 3.3.5 i R is equal to the output current of the Split-Pi converter which can be considered as input either to the load or to the central transmission network; in the case of the loads, the input voltage, v s , is coming from the central transmission network and at the output the load is considered as a resistor, R (Fig. 3.3.5). For the PV system, power profiles generated by the PV model in section 2.4.1.1 are taken into account as inputs to the Split-Pi converter. On the other hand, concerning the loads, we consider profiles obtained from several open source projects as in Department of Energy Office of Energy Efficiency & Renewable Energy (2017). Therefore, it remains to study the flat representation of the Split-Pi converter PH model, with v s the input voltage coming from source and i R the output current, as in (2.4.10)-(2.4.11) given below, which was presented also in PH form in section 2.4.2.1: Before continuing with the flat output calculation of the Split-Pi converter, we firstly consider the following: five states (3.3.31a-3.3.31e); p 1sc (t), p 2sc (t), q 1sc (t), q 2sc (t), q 3sc (t), four control inputs (3.3.31a-3.3.31e):

ṗ1sc (t) = q 1sc (t) C 1sc - q 2sc (t) C 2sc (1 -d 1sc (t)), (3.3.31a) ṗ2sc (t) = q 2sc (t) C 2sc (1 -d 2sc (t)) - q 3sc (t) C 3sc , (3.3.31b) q1sc (t) = v s R 1sc - q 1sc (t) C 2sc R 1sc - p 1sc (t) I 1sc , (3.3.31c) q2sc (t) = p 1sc (t) I 1sc (1 -d 1sc (t)) - p 2sc (t) I 2sc (1 -d 2sc (t)), (3.3.31d) q3sc (t) = p 2sc (t) I 2sc -i R (t). ( 3
d 1sc (t), d 2sc (t), v s (t), i R (t)
and the circuit's parameters and the values of the resistors are considered as constants,

1/C 1sc = α, 1/C 2sc = β, 1/C 3sc = γ, 1/I 1sc = η, 1/I 2sc = µ, 1/R 1sc = ν, 1/R = ζ.
Therefore, the equations below are deduced from the state-space representation (

    ω 1 ω 2 ω 3 ω 4     = P 1,[0] Z + [0]       dp 1sc dp 2sc dq 1sc dq 2sc dq 3sc       , (3.3.35)
where

P 1,[0] = ηp 1sc µp 2sc 0 βq 2sc 0 , Z + [0] =   0 1 0 0 0 0 0 0 1 0 0 0 0 0 -1   .
Therefore, the tangent flat outputs, ω j , j = 1, 2, 3, 4, of the Split-Pi converter are as follows:

ω 1 = ηp 1sc dp 1sc + µp 2sc dp 2sc + βq 2sc dq 2sc , (3.3.37a
)

ω 2 = dq 2sc , (3.3.37b
)

ω 3 = dp 2sc , (3.3.37c
)

ω 4 = -dq 3sc . (3.3.37d)
Consequently, after the integration of (3.3.37a)-(3.3.37d), we obtain the flat output set written below:

z 1 = η p 2 1sc 2 + µ p 2 2sc 2 + β q 2 2sc 2 , (3.3.38a) z 2 = q 3sc + q 1b , (3.3.38b) z 3 = q 2b , (3.3.38c) z 4 = q 2sc . (3.3.38d)
In the following , we introduce the states and inputs of the Split-Pi converter in function of the flat outputs (3.3.38a-3.3.38d) and their derivatives:

p 1sc = 2 η (z 1 -µ z 3 2 2 -β z 2 2
2

), (3.3.39a)

p 2sc = z 3 , (3.3.39b)
q 1sc = 1 α ( 1 2 1 2 η (z 1 -µ z 3 2 2 -β z 2 2 2 ) 2 η ( ż1 -µz 3 ż3 -βz 2 ż2 ) + βz 2 ż2 + µz 3 ż3 -γz 4 βz 2 η( 2 η (z 1 -µ z 3 2 2 -β z 2 2
2

)) , (3.3.39c)

q 2sc = z 2 , (3.3.39d) q 3sc = -z 4 , (3.3.39e
)

d 1sc = 1 - ż2 + µz 3 ż3 -γz 4 βz 2 η( 2 η (z 1 -µ z 3 2 2 -β z 2 2
2

)) , (3.3.39f) Furthermore, we consider the compact flat representation of the states and inputs in function of the flat outputs:

d 2sc = 1 - ż3 -γz 4 βz 2 , (3.3.39g) v s = 1 ν q1sc + αq 1sc + η ν p 1sc , (3.3
p 1sc = Φ 23 (z 1 , z 2 , z 3 ), (3.3.40a)
p 2sc = Φ 24 (z 3 ), (3.3.40b)
q 1sc = Φ 25 (z 1 , ż1 , z 2 , ż2 , z 3 , ż3 , z 4 ), (3.3.40c)
q 2sc = Φ 26 (z 2 ), (3.3.40d)
q 3sc = Φ 27 (z 4 ), (3.3.40e)
d 1sc = Φ 28 (z 1 , z 2 , ż2 , z 3 , ż3 , z 4 ), (3.3.40f) d 2sc = Φ 29 (z 2 , ż3 , z 4 ), (3.3.40g) v s = Φ 30 (z 1 , ż1 , z1 , z 2 , ż2 , z2 , z 3 , ż3 , z3 , z 4 , ż4 ), (3.3.40h) i R = Φ 31 (z 3 , ż4 ). ( 3

PV and loads flat representation through another set of flat outputs

This subsection proposes another set of flat outputs, taking into account the Proof 3.2, choosing four different states of the Split-Pi converter, since we have four control inputs:

z = p 1sc p 2sc q 2sc q 3sc . (3.3.41)
Consequently, through calculation, we obtain the states and the inputs in function of the aforementioned flat outputs:

p 1sc = z 1 , (3.3.42a)
p 2sc = z 2 , (3.3.42b) q 1sc = 1 α ( ż1 + βz 3 ż3 + µz 2 ż2 + γz 4 βz 3 ηz 1
), (3.3.42c)

q 2sc = z 3 , (3.3.42d) q 3sc = z 4 , (3.3
.42e)

d 1sc = 1 - ż3 + µz 2 ż2 + γz 4 βz 3 ηz 1 , (3.3.42f) d 2sc = 1 - ż2 + γz 4 βz 3 , (3.3.42g) v s = 1 ν q1sc + αq 1sc + η ν p 1sc , (3.3.42h) i R = µz 2 -ż4 . (3.3.42i)
Below, the corresponding compact flat representation of the states and inputs is introduced:

p 1sc = Φ 32 (z 1 ), (3.3.43a) p 2sc = Φ 33 (z 2 ), (3.3.43b) q 1sc = Φ 34 (z 1 , ż1 , z 2 , ż2 , z 3 , ż3 , z 4 ), (3.3.43c) q 2sc = Φ 35 (z 3 ), (3.3.43d) q 3sc = Φ 36 (z 4 ), (3.3.43e
)

d 1sc = Φ 37 (z 1 , z 2 , ż2 , z 3 , ż3 , z 4 ), (3.3.43f
) Following the previously depicted results for the ES system and the PV and loads system, we deduce that the less derivatives are interfering in the flat representation of the system, the less disturbances will emerge. Hence, in the upcoming section, we will concentrate on the concept of finding the proper set of flat outputs of a PH system considering the less possible derivatives. Furthermore, this information is useful for the B-spline parametrization, which will be used later in chapter 4, a method of parametrization which is strongly dependent on the order of derivatives appeared in the flat representation.

d 2sc = Φ 38 ( ż2 , z 3 , z 4 ), (3.3.43g) v s = Φ 39 (z 1 , ż1 , z1 , z 2 , ż2 , z2 , z 3 , ż3 , z3 , z 4 , ż4 ), ( 3

Port-Hamiltonian systems and differential flatness

Hereinafter, a connection between PH systems and differential flatness is presented. We propose a method for generating the possible sets of flat outputs for a particular class of PH systems.. The inverse dynamics of the PH systems will be described with illustrated examples considering a notion called bicausality within a Bond graph.

Bicausality and Bond graphs

As has been already mentioned in section 2.3.1, the causality property indicates the computational dependence among the ports characterized by pairs of effort e and flow f . Fig. 3.4.3 depicts the dependence on the efforts and flows between two ports in a causal Bond graph (for further details see also section 2.3.1). The red dashed lines illustrate the power flow and the dependence among the ports created after the causality assignment. Afterwards, another notion used in the Bond graph theory will be implemented. Bicausality is an extended version of causality which allows the interpretation of the inverse dynamics of a physical system [START_REF] Gawthrop | Bicausal bond graphs[END_REF][START_REF] Richard | Analysis of flatness using bond graphs and bicausality[END_REF]. In order to design a bicausal Bond graph, supplementary elements, such as the SS (source-sensor), should be taken into account. This elements enable the decoupling of ef pairs (Fig. 3.4.5) and allow the calculation of further quantitative variables in the Bond graph model. Within a bicausal Bond graph, we can analyze explicitly the inverse dynamics of a physical system concluding to the states or parameters estimation [START_REF] Ngwompo | Bond graph-based simulation of non-linear inverse systems using physical performance specifications[END_REF]). In the literature, a method exists called Sequential Causality Assignment Procedure for Inversion (SCAPI) [START_REF] Ngwompo | Bond graph-based simulation of non-linear inverse systems using physical performance specifications[END_REF], which describes the design process of the bicausal Bond graph. According to the SCAPI method, we follow the steps mentioned below: firstly, we define the minimal number of causal paths or the shortest causal path between the system input and output (Fig. 3.4.6); then, we replace each input and output element with the SS elements as in Fig. 3.4.6. The SS elements contain: the effort sources, S e , and the flow sources, S f , as inputs to the system, and the effort detectors, D e , and the flow detectors, D f , representing the outputs of the system. The sources elements have specific assignment as in Fig. 3 Note that for the junctions 0 and 1, the same rules are applied as in the case of the causal Bond graph (see also section 2.3.1).

f 6 = f 7 . (3.4.1f)
In the following, the differential causality of the storing and dissipative elements is defined:

f 5 = C de 5 dt = q, ( 3.4.2)
e 7 = I df 7 dt = ṗ, (3.4.3)
e R = f R R = (f 5 + f 6 )R = ( q + p I )R. (3.4.4)
where the integral causality can be taken also into account:

e 5 = 1 C t f 5 dt = q C , ( 3
p = If 7 = I( q + p I -q) = p, (3.4.7a)
q = Ce 5 = C ṗ, (3.4.7b) u = v s = ( q + p I )R + ṗ. (3.4.7c)
The state vector of the system is x = p q ∈ R 2×1 and the input vector is u = v s ∈ R 1×1 . Regarding the relations above, (3.4.7a)-(3.4.7c), through the bicausal Bond graph, a physical representation is obtained similar to the flat representation of the system after considering the flat outputs in function of the states (Definition 3.2.1). According to [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF], the number of flat outputs must be equal to the number of control inputs. The RLC system (3.4.7a-3.4.7c) has one control input, u = v s . Therefore, one flat output is necessary. From the bicausal Bond graph (Fig. 3.4.9), we retrieve the state p and obtain (3.4.7a), which shows that p cannot be written in function of the other variables of the RLC circuit, the state q and the input v s . Consequently, a straightforward choice is state p as a possible flat output for the system, z = p.

Therefore, the equations (3.4.7a)-(3.4.7c) are rewritten in function of the flat output as follows: Therefore, we concluded to a flat representation of the system through the bicausal Bond graph. Consequently, the bicausal Bond graph seems to be a suitable tool to analyze differential flatness, since both methods deal with the inverse dynamics of the system.

x = p q = z C ż , ( 3

Flat representation analysis of the ES system from bicausality

Hereinafter, through the ES electrical circuit, a method to derive the possible sets of flat outputs of an electrical circuit will be described. The electrical circuit of Fig. 3.3.1 is divided into two circuits: i) the KiBaM circuit with the resistance R 1b ; ii) the Split-Pi converter with the resistance R 1sc . The possible flat outputs of their PH representations will be generated through their bicausal Bond graphs comparing with the results of the previous sections. The idea, here, is based on the following proposition: Remark 3.4.1. (Link between bicausal Bond graphs and Dirac structures): If effort and flow variables of the bicausal Bond graph are related with equations as in (2.5.5) of Definition 2.5.1, then the system is written in kernel representation of a DS. Hence, the overall dynamics represents a PH system. As a consequence, in the case of electrical circuits, the flat representation of their PH models can be analyzed through their bicausal Bond graphs as already presented for the RLC circuit in section 3.4.1. Indeed taking into account the bicausal Bond graphs (Fig.

+ F f I1 f I2 . . . f In I f C1 f C2 . . . f Cn C f R1 f R2 . . . f Rn R y ,
where n C is the number of capacitors C, n I is the number of inductors I (with the capacitors and inductors being the storing elements) and n R is the number of the dissipative elements. The variable u is the input and the variable y is the output of the system. Therefore, if the kernel representation exists, it composes a DS, hence a PH system, according to Definition 2.3.3, where a PH system is expressed as noted below:

(-ẋ, ∇H(x), f R , e R , f E , e E ) ∈ D. ( 3
q 2 C + 1 2
p 2 I (see also section 2.4.1). Therefore, the flat representation of the corresponding PH system can be analyzed as in section 3.4.1 for the RLC circuit.

KiBaM battery

Firstly, we concentrate on the battery (see Fig. 3.4.10) of which the PH state-space representation, presented also in section 2.4.1.2 in PH form, is: where the state vector is x b = q 1b q 2b ∈ R 2×1 and the input vector u b = v sc ∈ R 1×1 , with v sc the voltage coming from the Split-Pi converter. According to flatness theory, since we have one input, we have to find one flat output. Developing the inverse dynamics of the battery's system through its bicausal Bond graph, we will deduce the possible flat outputs of the battery. Fig. 3.4.11 presents the bicausal Bond graph of the battery, according to the SCAPI method previously described. The graph is composed by two storing elements, the capacitors C 1b and C 2b with their differential causalities, one effort source, S e , which is the input, and one flow detector, D f , which is the output. Therefore, the state and input relations of the inverse system are written below:

q1b = - 1 R 1b + 1 R 2b q 1b C 1b + 1 R 2b q 2b C 2b + 1 R 1b v s , (3.4.12a) q2b = 1 R 2b q 1b C 1b - 1 R 2b q 2b C 2b . ( 3
q 1b = C 1b f 5 = q2b R 2b C 1b + C 1b C 2b q 2b , ( 3
.4.13a) As a next step, we confirm if the system obtained by the bicausal Bond graph is a DS, hence a PH system. To demonstrate this, the kernel representation of the system is described which, according to Definition 2.5.1 (Ee + F f = 0), if present (verifying the conditions EF + F E = 0 and rank(F + E) = dimF ), then the system composes a DS. Therefore, the kernel representation of the bicausal Bond graph in Fig. 3.4.11 is presented below:

q 2b = C 2b f 9 = v s C 2b -q1b R 1b C 2b -q2b R 1b C 2b -q2b R 2b C 2b , (3.4.13b) v s = q1b R 1b + q2b R 1b + q2b R 2b + q 2b C 2b . ( 3
      1 -1 0 1 0 -1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0       E       v q 1b v q 2b v R 1b v R 2b v sc       e +       0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 -1 0 1 0 1 0       F       i q 1b i q 2b i R 1b i R 2b i s       f =       0 0 0 0 0       , ( 3
q 2 C + 1 2 p 2 I .
Consequently, all the possible relations among the efforts and flows, derived from the 0 and 1 junctions (Fig. 3.4.11), are presented below:

        q 1b q 1b q 2b q 2b v sc v sc         =            0 0 C 1b C 2b R 2b C 1b 0 0 -R 1b C 1b 0 -R 1b C 1b C 1b 0 0 1 0 0 0 0 0 0 0 0 R 1b 1 C 2b R 1b + R 2b 0 0 0 0 0 0                  q 1b q1b q 2b q2b v sc       . (3.4.15)
From the latter matrix representation, in the left part the states and the inputs are illustrated, which, concerning flatness theory, can be written in function of the flat outputs and their derivatives. Since we have one control input, v sc , we expect only one flat output [START_REF] Levine | Analysis and control of nonlinear systems: A flatness-based approach[END_REF]. Taking into account Proof 3.2, where the flat output can be written in function of the states, we assume that the only possible flat output for the battery is z = q 2b from (3.4.15). This conclusion comes from the fact that it is the only state not written in function of the inputs, the other states or their derivatives. In order to look at the result more carefully, we separate the matrices extracting the circuit parameters, C 1b , C 2b , and the dissipation, R 1b , R 2b , as follows:

           q 1b q 1b C 1b q 2b C 2b q 2b C 2b u s u s            =                     0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 1         V +         0 0 0 R 2b 0 0 -R 1b 0 -R 1b 0 0 0 0 0 0 0 0 0 0 0 0 R 1b 0 R 1b + R 2b 0 0 0 0 0 0         W                     q 1b C 1b q1b q 2b C 2b q2b v sc         . ( 3

.4.16)

Therefore, from matrix V we can derive the states and control inputs which are structurally required for the flat representation. The state q 2b is included, in function of which q 1b and v sc can be derived (third column of matrix V ). In W matrix, the dissipative elements are presented and are involved only with the derivatives of the states, which are not taken into consideration for the flat outputs selection. The flat representation of the battery is written below from (3.4.13a)-(3.4.13c):

x b = q 1b q 2b =   R 2b C 1b ż + C 1b C 2b z z   , (3.4.17a) u = v sc = R 2b R 1b C 1b z + (R 1b C 1b C 2b + R 1b + R 2b ) ż + 1 C 2b z. (3.4.17b)
The compact flat representation, as in Definition 3.2.1, of the KiBaM battery circuit deriving from (3.4.17a) and (3.4.17b) is:

q 1b = Φ 44 (z, ż), ( 3
.4.18a)

q 2b = Φ 45 (z), (3.4.18b) v sc = Φ 46 (z, ż, z). (3.4.18c)
This approach can be extended for the nonlinear cases of constitutive equations.

Split-Pi converter

Afterwards, following the same steps as for the KiBaM battery, the flat representation of the Split-Pi converter will be analyzed through its bicausal Bond graph. At first, its PH state-space representation is presented (see Fig. 3.4.10) below, introduced also before:

ṗ1sc = q 1sc C 1sc - q 2sc C 2sc (1 -d 1sc ), (3.4.19a) ṗ2sc = q 2sc C 2sc (1 -d 2sc ) - q 3sc C 3sc , (3.4.19b) q1sc = v s R 1sc - q 1sc C 1sc R 1sc - p 1sc I 1sc , (3.4.19c) q2sc = p 1sc I 1sc (1 -d 1sc ) - p 2sc I 2sc (1 -d 2sc ), (3.4.19d) q3sc = p 2sc I 2sc -i R . (3.4.19e)
where the state vector is x sc = p 1sc p 2sc q 1sc q 2sc q 3sc ∈ R 5×1 and the input vector

u sc = v s i b ∈ R 2×1
, where v s is the voltage coming from the source and i b is the current of the resistor R. Furthermore, the converter's duty cycles, d 1sc , d 2sc are considered also as control inputs to the system (see also Fig. 3.3.1).

Since we have four control inputs in total, we have to find four flat outputs. Developing the inverse dynamics of the system through its bicausal Bond graph, the four possible states considered as flat outputs will be deduced.

p 1sc = I 1sc u s R 1sc - q 1sc C 1sc R 1sc -q1sc , (3.4.20a) p 2sc = I 2sc [ q3sc + i R ] , (3.4.20b) q 1sc = C 1sc ṗ1sc + q 2sc C 2sc (1 -d 1sc ) , (3.4.20c) q 2sc = q 2sc , (3.4.20d) q 3sc = C 3sc q 2sc C 2sc (1 -d 2sc ) -ṗ2sc . (3.4.20e)
Moreover, the control inputs are deduced also in function of the states and their derivatives from the bicausal Bond graph as follows:

v s = p 1sc I 1sc + q 1sc C 1sc R 1sc + q1sc R 1sc , (3.4.21a) i R = q 3sc C 3sc R , (3.4.21b) d 1sc = C 2sc q 2sc ṗ1sc - q 1sc C 1sc + q 2sc C 2sc , (3.4.21c) d 2sc = C 2sc q 2sc -ṗ2sc + q 2sc C 2sc - q 3sc C 3sc . (3.4.21d)
Next, the kernel representation of the bicausal Bond graph of the Split-Pi converter is presented (Definition 2.5.1), which proves that it composes a DS:

            -1 0 -1 1 -d 1sc 0 0 0 0 0 1 0 1 -d 2sc -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 0 1 0 0 0 0 -1 0 -1 0 0 0 0 0 0 0 0 0             E             e p1sc e p2sc e q1sc e q2sc e q3sc e R1sc e R v s             e + (3.4.22a) +             0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 1 0 0 0 0 1 1 -d 1sc -(1 -d 2sc ) 0 1 0 0 0 0 0 -1 0 0 -1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1             F             f p1sc f p2sc f q1sc f q2sc f q3sc f R1sc f R i s             f =             0 0 0 0 0 0 0 0             , Consequently,
q 2 C + 1 2 p 2 I
). Subsequently, all the relations among the flow and the effort variables from the 0 and 1 junctions (Fig. 3.4.13) are introduced below:

p 1sc p 1sc p 2sc p 2sc q 1sc q 1sc q 2sc q 3sc q 3sc q 3sc v s i R = =                             1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - I 1sc C 1sc R 1sc -I 1sc 0 0 0 I 1sc R 1sc 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 I 2sc 0 I 2sc 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 C 1sc 0 0 0 0 C 1sc C 2sc 0 0 0 0 -C 1sc 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 -C 3sc 0 0 C 3sc C 2sc 0 0 0 0 0 -C 3sc 0 0 0 0 0 0 0 0 0 0 C 3sc R 0 0 R 1sc I 1sc 0 0 0 1 C 1sc R 1sc 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 C 1sc R 0 0 0 0 0                             • • p 1sc ṗ1sc p 2sc ṗ2sc q 1sc q1sc q 2sc q 3sc q3sc v s i R d 1sc d 2sc . (3.4.23)
At this point, four additional relations are obtained for the states p 1sc , p 2sc and q 2sc which include the duty cycles of the converter, d 1sc and d 2sc :

p 1sc I 1sc d 1sc = p 1sc I 1sc -q2sc - p 2sc I 2sc (1 -d 2sc ), (3.4.24a)
p 2sc I 2sc d 2sc = p 2sc I 2sc -q2sc - p 1sc I 1sc (1 -d 1sc ), (3.4.24b)
q 2sc C 2sc d 1sc = q 2sc C 2sc - q 1sc C 1sc + ṗ1sc , (3.4.24c)
q 2sc C 2sc d 2sc = q 2sc C 2sc - q 3sc C 3sc -ṗ2sc . (3.4.24d)
These four equations will be considered in the analysis of the flat outputs under restrictions. This is because of the fact that during the operation of the converter the states p 1sc , p 2sc , q 2sc are different from 0. Contrariwise, the duty cycles, d 1sc and d 2sc , can be equal to 0. Therefore, the duty cycles cannot be placed in the denominator's position and cannot be included in the flat outputs selection. Afterwards, the circuit's parameters and the dissipative elements are divided into two matrices as in (3.4.16), V and W and (3.4.23) becomes:

p 1sc p 1sc I 1sc p 2sc p 2sc I 2sc q 1sc q 1sc C 1sc q 2sc q 3sc q 3sc C 3sc q 3sc C 3sc v s i R = = (V + W ) p 1sc ṗ1sc p 2sc ṗ2sc q 1sc q1sc q 2sc q 3sc q3sc v s i R d 1sc d 2sc , (3.4.25)
where

V =                    
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 -1 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

                    (3.4.26)
and

W =                        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 1 R 1sc 0 0 0 0 1 R 1sc
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 R 0 0 R 1sc 0 0 0 R 1sc 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 R 0 0 0 0 0                        . (3.4.27)
From (3.4.23), we assume that the set of flat outputs for the Split-Pi converter system contains the states p 1sc , p 2sc , q 1sc , q 2sc , q 3sc . Concerning the state q 2sc , it must be regarded as one of the flat outputs (seventh row of matrix V ). Since it cannot be expressed in function of the other states and inputs, it is necessary for the flat representation. Additionally, because of (3.4.24a)-(3.4.24d), it is indispensable to consider states p 1sc , p 2sc , q 2sc as part of the flat outputs since they are linked to the switching activity of the converter. Moreover, the d 1sc and d 2sc are given only from (3.4.24a)-3.4.24d and they cannot be taken into account as flat outputs for the system. The control inputs, v s and i R , similarly cannot be considered, since they are written in function of the states and their derivatives, but not all the states can be written in function of v s and i R (eleventh and twelfth column of matrices V and W ).

Thereupon, referring to the noise that the derivatives of the states create in the results, we need to select the states appropriate from (3.4.25), (3.4.26), (3.4.27) and (3.4.23) (less derivatives create less deviations among the reference and the actual values), so that the less possible derivatives appear in the flat representation. If we look into (3.4.23), the state q 1sc is in function of more than one derivatives, while the state q 3sc can be replaced from an equation with no derivatives included. However, neither can be excluded and we conclude to the following possible flat outputs sets: z 1 = p 1sc p 2sc q 2sc q 1sc , (3.4.28) for the first flat output z 1 = p 1sc p 2sc q 2sc q 1sc , the flat representation of the states and the inputs is introduced: representation of the Split-Pi converter. Although, in DC microgrids, the voltage regulation of the central transmission network is one of the principal factors to investigate. Here, the input voltage, v s of the Split-Pi converter is considered in the optimization problem chapter 4, where the flat representation is combined with the B-spline parametrization, strongly dependent on the number of derivatives included. Consequently, to simplify the calculations at this case (see also (3.4.21a)), the first set of flat outputs (3.4.28) would be the appropriate choice.

z 2 = p 1sc p 2sc q 2sc q 3sc . ( 3
p 1sc = z 1 , (3.4.30a) p 2sc = z 2 , (3.4.30b) q 1sc = z 4 , (3.4.30c) q 2sc = z 3 , (3.4.30d) q 3sc = 1 γ βz 3 ż3 µz 2 + ηα β z 1 z 4 z 3 + ηz 1 ż1 -ż2 , (3.4.30e) d 1sc = 1 - µ ż3 z 2 η 2 z 2 1 γ β q 3sc z 3 -ż2 , (3.4.30f) d 2sc = 1 - ż3 z 2 µ -ηz 1 α β z 4 z 3 + ż1 , (3.4.30g) v s = 1 ν ηz 1 + ναz 4 + ż4 , (3.4.30h) i R = µz 2 -q3sc . ( 3

Conclusion

In this chapter, the differential flatness method was introduced and implemented for PH systems. The flat representation was investigated for the various components (ES, PV, loads) of the reduced DC microgrid presented in PH form in chapter 2. The flat outputs were, primarily, generated through an algorithm proposed in [START_REF] Franke | On the computation of flat outputs for nonlinear control systems[END_REF]) which provides the flat representation of the tangent system. However, since multiple sets of flat outputs can exist for each physical system, to compare and be able to decide what is the best choice, another flat output was found according to Proof 3.2, where the flat outputs can be written in function of the states of the system. The sets of flat outputs were validated through simulations, where the noise appeared in the flat system, created by the numerical differentiation, affects the signals response. The conclusion was that the more derivatives exist in the flat representation of the states and control inputs, the more noise-affected the results are. Afterwards, the concept of bicausality in the Bond graphs was described which considers the inverse dynamics of physical systems. Through several examples, it was shown that bicausality is an appropriate tool to study flat representation. This was further extended in PH systems proving through the kernel representation that the bicausal Bond graphs compose a Dirac structure. Afterwards taking into account all the possible relations among the pairs of efforts and flows, a primary procedure was described which generated in a straightforward way the possible flat outputs of the illustrated examples.

In the succeeding chapter, the flat representation of the system will be employed in the energy management problem considered for the DC microgrid and will be combined with B-spline parametrization for continuoustime constraint validation.

Chapter 4

Multilevel supervision of a DC microgrid

Introduction

In the previous chapters, the proposed DC microgrid system was built stepwise, describing its detailed dynamical model using PH modeling. Next, through differential flatness and bicausality properties, the inverse dynamics of the system components (PV, ES, loads) were provided. In this chapter, the global and local supervision of the system will be investigated. To select a suitable control approach, various issues need to be taken into consideration, as for example the multiple timescales or the islanding and grid-connecting modes. As already mentioned in chapter 1, a variety of control approaches was developed according to the different requirements appearing in the energy management problem.

In general, the energy management problem of a microgrid is considered as a constrained optimizationbased control problem not straightforward to solve. Considering that fast (converters) and slow (DERs) dynamics components which are composing the system, a control system in many levels needs to be implemented. The higher level with longer planning intervals and global supervision tasks will account for the lower level with smaller time intervals and local supervision. The most suitable control architecture to deal with these problems is, consequently, the hierarchical control. Its design framework is organized as a 'tree', as already introduced in chapter 1 in Fig. 1.1.1, from the higher level up to the lower level. In the literature, various control strategies were proposed to date, where the hierarchy presides (Bouzid et 2018), they present a combination of MPC with machine learning, using MPC for the optimal profile generation and machine learning for system adjustment to the consumer's demand; the secondary control (middle level) enforces the amelioration of power quality. It targets the internal processes of the system under voltage and frequency disturbances. Both centralized [START_REF] Shafiee | Distributed secondary control for islanded microgrids-a novel approach[END_REF] and distributed [START_REF] Lou | Distributed mpc-based secondary voltage control scheme for autonomous droop-controlled microgrids[END_REF] secondary controllers have been investigated and various control methods have been proposed. In [START_REF] Shafiee | Distributed secondary control for islanded microgrids-a novel approach[END_REF], a phase controller is used which regulates the phase angle of the distributed generators instead of the frequency to decrease frequency and amplitude discrepancies. In [START_REF] Wang | An improved distributed secondary control method for dc microgrids with enhanced dynamic current sharing performance[END_REF], the authors aim at reducing the voltage deviations caused by droop control method in a distributed control framework for accuracy in current sharing. [START_REF] Baghaee | Performance improvement of multi-der microgrid for small-and large-signal disturbances and nonlinear loads: novel complementary control loop and fuzzy controller in a hierarchical droop-based control scheme[END_REF] improves the performance of a microgrid by controlling the reactive power under disturbances caused by power outages, short circuits and the like. In Simpson-Porco et al. ( 2015), frequency and voltage controllers are designed, enabling the possibility for the system to achieve either frequency or voltage regulation while maintaining reactive power sharing. Predictive controllers are used such as MPC [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]. For instance, [START_REF] Ahumada | Secondary control strategies for frequency restoration in islanded microgrids with consideration of communication delays[END_REF] implements MPC in combination with a PI controller in order to decrease the discrepancies between the nominal and the actual frequency, to increase the stability of the system and to counteract communication delays; the primary control (low level) includes a localized supervision of the power distribution and the voltage/current adjustment among the DERs and converters. This control level defines the proper operation of the converters in order to ensure stability. The converters contain an internal switching activity that obeys to an external power loop based on a management strategy [START_REF] Hadjeras | Hybrid adaptive control of the boost converter[END_REF][START_REF] Wang | Hybrid energy storage system using bidirectional single-inductor multiple-port converter with model predictive control in dc microgrids[END_REF]. Such strategies are for instance the maximum power point tracking [START_REF] Liu | A review of maximum power point tracking methods of pv power system at uniform and partial shading[END_REF], which contains controllers that extract the maximum available power from energy sources, like wind turbines or PVs. Another approach is the master/slave control [START_REF] Mazumder | Master-slave current-sharing control of a parallel dc-dc converter system over an rf communication interface[END_REF], where a device is selected as the master controller over one or more devices that act as slaves. Additionally, there is also droop control [START_REF] Guerrero | Hierarchical control of droop-controlled ac and dc microgrids-a general approach toward standardization[END_REF], that allows balance of load sharing during parallel generators operation. Moreover, inner voltage and current control loops are further considered to study the stability of the system (Gavagsaz-Ghoachani et al., 2016).

Another important aspect which influences the operation of microgrids is the power loss. Optimal power flow accounts for the effective and reliable functioning of the system with either minimizing the electricity cost or the energy dissipation. Power losses in the DC-bus network can significantly affect the power quality during transmission. Hence, their mitigation constitutes a very important factor for the improvement of the power transmission. Different approaches have been proposed recently for the power loss reduction, concentrating either on the components connected to the microgrid or on the central transmission network [START_REF] Gamarra | Computational optimization techniques applied to microgrids planning: A review[END_REF]. Some works focus on topological issues or optimal scheduling of the energy storage (ES) and the renewable sources, as in [START_REF] Iovine | Nonlinear control of a dc microgrid for the integration of photovoltaic panels[END_REF] and [START_REF] Wei | On optimally reducing power loss in microgrids with power storage devices[END_REF]. Others concentrate on the existing power losses in the transmission lines, which are inherent to the electrical network and they cannot be prevented or eliminated. In the literature, researchers have proposed several methods for optimal power distribution including power loss minimization. In [START_REF] Nahata | Hierarchical control in islanded dc microgrids with flexible structures[END_REF], the authors propose a three-layer hierarchical control approach to solve the energy management problem in islanded microgrids. More specifically, at the higher level, taking into account topological and stability issues, optimal power profiles are generated, within a MPC framework. At the middle level, using the upper level references, a voltage regulation problem is solved taking into account the power losses of the transmission network and the converters. Furthermore, a two-level hierarchical control based on plug and play method has been proposed in [START_REF] Vazquez | A fully decentralized adaptive droop optimization strategy for power loss minimization in microgrids with pv-bess[END_REF] where at a low level load sharing is achieved through droop control and at a higher level adaptive droop controllers are developed for power losses minimization. [START_REF] Sanseverino | Optimal power flow in islanded microgrids using a simple distributed algorithm[END_REF] introduce an optimal power flow algorithm for islanded microgrids to cope with central power dissipation under a centralized three-layer supervision controller.

In the methods and techniques previously described, the dynamical models are mostly considered as sets of differential equations [START_REF] Parisio | Stochastic model predictive control for economic/environmental operation management of microgrids: An experimental case study[END_REF]. Consequently, the models are described as sets of constitutive equations which do not explicitly represent the structure of the power-preserving functions or maintain the energy conservation within the system. Next, concerning the generation of the optimal profiles, constrained optimization-based control methods are often employed such as MPC [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF], a popular method to generate on-line optimal profiles for discrete-time systems. Depending on the type of MPC, the cost function penalizes the cost, the dissipation or the error among the actual and reference signal profiles. However, the microgrids are convoluted networks where numerous factors need to be examined and considered at the same time, such as power optimization, cost minimization, stability, robustness and the like. Consequently, no definitive microgrid models and control methods exist and many aspects remain still open for investigation. In this work, a different approach is developed in a hierarchical control framework. The models already introduced in chapter 2 will be considered in PH formulation, which always respects the power conservation among the physical elements and describes the dynamics of the system explicitly. Furthermore, differential flatness for continuous-time optimal profile generation will be used, combined with B-spline approximation for continuous-time constraint validation.

This chapter presents the power balancing problem and aims at reducing the electricity cost and the power dissipation exploiting the energy storage and the meshed topology of the system. The contributions of this chapter are summarized below: the use of differential flatness and B-splines parametrization at the high level to obtain optimal reference profiles for the UG, ES power and for the battery's voltage and current. In contrast to the previously referenced publications, which concentrate more on obtaining optimal profiles for the source power generation, in here a further analysis with optimal profiles in continuous time is provided not only for the power but also for the current and voltage of the battery, which will be used at the middle level as reference. Through B-spline approximation of the flat representations, the continuous-time constraint validation will be ensured; the multi-level control structure which provides at each level (high, middle and low) optimal profiles to be followed by the lower level. For instance, power balancing requirements at the high level lead to an optimal battery usage profile, to be tracked at the middle level; the middle level provides voltage and current references for the battery which are tracked at the low level by an explicit switching law for the DC/DC converters. In all cases bottom to top information in the power flow is exploited for a reliable profile generation (e.g., tracking errors are accounted for); the validation of the proposed hierarchical control through extensive simulations based on realistic load, renewable power and electricity price profiles. The behavior of the scheme is analyzed at the high, middle and low levels over the meshed DC microgrid benchmark previously introduced. This chapter is organized as follows. Firstly, the basic tools used to build the hierarchical control problem are introduced. Then, the energy management problem will be presented defining, primarily, the objectives, the reference profiles and the set of constraints. Afterwards, the multi-layer supervision problem is analyzed explicitly. Finally, the simulation results and conclusions are presented.

Constrained optimization-based problem formulation

Optimal control theory deals with systems which can be controlled under a specific optimal criterion with given inputs, boundary conditions and a set of constraints for the states. It's a method that defines the evolution of a dynamical system generating state trajectories over a period of time. This section recalls the basic tools important to construct the multi-scale supervision problem for the system presented in chapters chapter 2 and chapter 3. The notions of optimal control theory, B-spline parametrization and MPC tracking are presented [START_REF] Rawlings | Model predictive control: Theory and design[END_REF][START_REF] Suryawan | Constrained Trajectory Generation and Fault Tolerant Control Based on Differential Flatness and B-splines[END_REF].

Constrained optimization-based control

A constrained optimization problem for a control system as in (3.2.1) contains an objective function and a set of constraints as follows [START_REF] Bertsekas | Constrained optimization and Lagrange multiplier methods[END_REF]. The general constrained minimization problem is written as follows:

min u J (x(t), u(t)) (4.2.1) subject to ẋ(t) =f (x(t), u(t)), x(t 0 ) = x 0 ∀t ∈ [t 0 , t f ], (4.2.2)
c eq (x(t), u(t)) = 0, eq = 1, . . . , N ceq , ∀t ∈ [t 0 , t f ], (4.2.3)

c neq (x(t), u(t)) ≤ 0, neq = 1, . . . , N neq , ∀t ∈ [t 0 , t f ], (4.2.4) (t f , x(t f )) ∈ S f . (4.2.5)
The objective function (otherwise called cost function) J is defined as:

J (x(t), u(t)) = φ 1 (x t f ) + t f t0 φ 2 (x(t), u(t))dt, (4.2.6)
with [t 0 , t f ] the time interval of interest, φ 1 (x t f ) and φ 2 (x(t), u(t)) are given functions, the terminal and the running cost, respectively, and x t f is the final state vector. The cost function is chosen based on different objectives, as the electricity cost minimization or the power dissipation. In addition to this, the states and control variables obey to a set of equality and inequality constraints referred below:

c eq (x(t), u(t)) = 0, eq = 1, . . . , N ceq , ∀t ∈ [t 0 , t f ], c neq (x(t), u(t)) ≤ 0, neq = 1, . . . , N neq , ∀t ∈ [t 0 , t f ], (4.2.7)
with c eq the equality constraints, c neq the inequality constraints and N eq , N neq the number of equality and inequality constraints correspondingly. These relations arise from the system requirements, parameters and some general physical restrictions. Furthermore, there is another set of constraints, S t f , which defines the endpoints of the final time t f and the final state x t f which may vary depending on the optimization problem. Consequently, the final time and final state must additionally fulfill the constraint below:

(t f , x(t f )) ∈ S f . (4.2.8)
In general, a constrained optimization-based control problem is not straightforward to solve, especially when there are nonlinearities in the objective functions and the constraints. Usually the solution to these problems is given through methods which require system discretization. Contrarily to the classical methods, in here the optimization problem (4.2.1) is maintained in continuous time, applying the notions of differential flatness from chapter 3 and B-spline parametrization, described below.

B-spline parametrization

In this section, the B-spline parametrization and some of its basic properties will be presented. The flat output will be projected over a finite set of basis functions, called B-splines, in function of which the cost and constraints will be written later in this chapter. The idea of B-spline parametrization is based on the construction of a curve defined by a number of points. These points are connected successively representing the polynomials which characterize and form the B-splines (see also Fig

i,1 (t) = 1, τ i ≤ t < τ i+1 , 0, otherwise, b i,d (t) = t -τ i τ i+d-1 -τ i b i,d-1 (t) + τ i+d -t τ i+d -τ i+1 b i+1,d (t). (4.2.9)
The b-splines are defined by a knot-vector T ∈ R N +d which is a set of non-decreasing time instants equal to

T = {τ 0 ≤ τ 1 ≤ ... ≤ τ κ }, where τ κ = N + d + 1: T = {τ 0 , τ 1 , . . . , τ κ-1 , τ κ } = {t 0 , . . . , t 0 d , τ d , . . . , τ κ-d-1 , t f , . . . , t f d }, (4.2

.10)

where τ 0 = t 0 , τ κ = t f .

z(t) = N i=1 p i • b i,d (t) = PB d (t), (4.2

.11)

where p i is the i th control point, to which it corresponds the matrix P of N control points equal to

P = p 1 p 2 ... p N ∈ R N ×1 . Moreover, B d (t) = b 1,d (t) b 2,d (t) ... b N,d (t) ∈ R N ×d is the B-spline vector.
B-spline parametrization seems appropriate due to its numerous properties such as convexity, continuity, ease on derivatives computation and so on. The B-splines are used because they can easily enforce continuity across the control points. This accounts for continuous-time constraint validation considering the control points as decision variables [START_REF] Stoican | Constrained trajectory generation for uav systems using a b-spline parametrization[END_REF]. Furthermore, they facilitate the computation of the derivatives existing from the flat representation. Below, some of the most significant properties of B-splines are enumerated (Stoican et P2. The r-order derivatives of d -order B-splines can be expressed as linear combinations of dr order B-splines which, in turn, can be expressed as d -order B-splines over each knot sub-interval:

B (r) d (t) = M d ,d -r B d -r (t) = M d ,d -r S κ,d -r,d B d (t), ∀t ∈ [τ κ , τ κ+1 ), (4.2.12)
with S κ,d -r,d denoting the translation matrix from higher to lower degree basis functions and M d ,d -r the matrix performing the linear combinations of the lower-degree basis functions.

P3. The sum of the B-splines basis functions b i,d (t) is equal to 1 for the time interval [τ 0 , τ κ ] and the B-splines themselves are either equal or greater than 0 (Stoican et al., 2017; Suryawan, 2012).

P4. Every B-spline function depends on its order r and its degree d . The degree of B-splines depends on the derivative order where the continuity is ensured. Each B-spline may be locally treated without influencing the rest of the curve.

P5. From the third property, we have that for the interval [τ κ , τ κ+1 ] the sum of B-splines basis functions b i,d (t) is equal to 1. Therefore, from (4.2.11)and (4.2.9), z(t) can be written also as:

z(t) = i:t∈[τκ,τκ+1] i j=i-d +1 p j • b i,d (t) (4.2.13)
and exists in the following unity of all convex hulls formulated by the control points (Stoican et al., 2017; Suryawan, 2012):

z(t) ∈ n i=d -1 conv{p i-d +1 , . . . , p i }, (4.2.14)
where n is the total number of convex hulls.

Model Predictive Control

An important optimization-based control method on which this thesis builds its results is Model Predictive Control (MPC). We will consider two types of MPC, tracking MPC and economic MPC. The tracking MPC will be taken into account in the middle level for reference tracking and the economic MPC will be considered as a comparison tool for generating optimal profiles at the high level.

MPC is an on-line optimal control technique which optimizes a cost function while satisfying a set of constraints. The tracking MPC significantly reduces the deviations among the real and reference profiles, while the economic MPC penalizes a general cost depending on the objective, such as electricity cost or power loss mitigation [START_REF] Ellis | A tutorial review of economic model predictive control methods[END_REF][START_REF] Rawlings | Model predictive control: Theory and design[END_REF]. In general, MPC algorithms take into account discretized dynamical systems, the current state and a cost function over a receding prediction horizon in order to compute the future evolution of the system.

The MPC controller solves an open-loop optimization problem at each time step k over a prediction horizon (see also the proof of concept illustration given in Fig. 4.2.2). Next, it applies the first value of the obtained control sequence to the system dynamics, updates the system's states and repeats the procedure at the next step. Receding horizon control (RHC) was one of the first attempts done by researchers [START_REF] Kamel | Model predictive control for trajectory tracking of unmanned aerial vehicles using robot operating system[END_REF][START_REF] Mayne | Receding horizon control of nonlinear systems[END_REF][START_REF] Rawlings | The stability of constrained receding horizon control[END_REF] to approach more closely problems with fixed horizon and allowed the computation of a complete control sequence by sending only the first step to the system in order to recompute and restart the process. Consider the following nonlinear discrete-time system:

x(k + 1) = f d (x(k), ũ(k)), x(t) = x, k ≥ t, (4.2.15)
with the state vector x(k) ∈ R n and the control input u(k) ∈ R m . For the considered optimization problem, the cost function to minimize is defined as follows (Allgöwer and Zheng, 2012):

J M P C = t+Np-1 k=t (x(k) Q x x(k) + ũ(k) R ũ ũ(k)) + Ṽf (x(t + N )), subject to : x(k + 1) = f d (x(k), ũ(k)), x(k) ∈ X , ũ(k) ∈ Ũ, x(t + N p ) ∈ Xf (terminal constraint),
where Q x and R ũ are matrices of appropriate dimensions, X is the set of state constraints and U the set of input constraints. Furthermore, the term Ṽf (x(t + N )) plays an important role for the system's stability, since it forces the states to always exist in a particular set.

Control objectives, reference profiles and general constraints

This subsection presents the ingredients necessary for the formulation of the flatness-based optimization control problem of the reduced meshed DC microgrid introduced in Fig. 2.5.1.

Control objectives

The energy management problem for the meshed DC microgird system (Fig. 2.5.1) addresses the following control objectives: energy management, load balancing in the common DC-bus among the nodes; cost minimization exploiting the PV power generation and the energy storing devices; reduction of the power losses during distribution within the central transmission network; satisfaction of a set of constraints for several variables of the system, such as current, voltage, charge etc.

continuity in the operation of the system in case of a line under fault.

Reference profiles

Fig. 4.2.3 presents the reference profiles of the electricity price, the PV panel and the loads which will be used later in the simulations:

a

Consumers demand reference profile for commercial use

Constraints

In this section, we will describe the necessary constraints and restrictions for the system as previously described in section 2.6.2. As already emphasized, the ES system is of high importance and its long lifetime and maintenance are significant factors. Because of the fact that the ES stores and conserves energy, it ensures the safe and continuous operation of the DC microgrid when other resources (PV and UG) are unable to generate power. Priority on the battery's optimal operation will be given and the dynamics of the ES system will be analytically described in the optimization problem. Therefore, setting constraints for the battery's proper operation is indispensable. Firstly, since the battery can store a finite amount of energy, boundaries concerning the charge must be considered. The charge of the battery has to be between a minimum q min b and a maximum q max b value (both states, q 1b (t) and q 2b (t) as in (2.4.8), have positive values and q max b = q max 1b + q max 2b . Additionally, q min b and q max b depends on the depth of discharge (DOD), which describes the battery's discharging possibilities.

12.7V

According to the characteristics of the AGM 12-165 lead acid batteries (Victron Energy, 2015), they are fully charged at 92%. Therefore, the DOD is equal to 8%. The maximum capacity of the battery is equal to 165Ah. By respecting the DOD of the battery, the battery will keep its rated capacity for a limited number of charge/discharge cycles, which are shown in Table 4.2.1. Furthermore, another important factor is the charging and discharging rate. A quick discharge can reduce the amount of generated power and can influence the battery's capacity. This is because of the limited time remaining for the battery's chemical elements to move to their necessary positions, which leads to the decrease of the available power inside the battery. In the contrary, a slow discharge rate contributes to the amelioration of the capacity and the generated power. The rate of charge is defined by the current of the battery which is also bounded between a maximum, i min b , and a minimum, i max b , value. Similar constraints exist for the voltage of the battery.

All the aforementioned constraints for the proper operation of the battery are expressed below (see also Table 4.2.2:

q min b ≤ q b (t) ≤ q max b , (4.2.17a)
Charging mode : (4.2.17b)

v min b,charging ≤ v b,charging (t) ≤ v max b,charging , (4.2.17c)
i min b,charging ≤ i b,charging (t) ≤ i max b,charging , (4.2.17d)
Dischaging mode : (4.2.17e)

v min b,dicharging ≤ v b,dicharging (t) ≤ v max b,dicharging , (4.2.17f)
I b,dicharging (t) 15A. (4.2.17g)
In addition, boundaries are defined for the voltage in the central transmission network (Fig. 2.6.2), which maintain the voltage on the connecting nodes close to 400 V , as already presented in chapter 2:

v min,h DC ≤ v ug (t), v pv (t), v es (t), v loads (t) ≤ v max,h DC , (4.2.18)
where v min,h DC is equal to 380 V and v max,h DC is equal to 430 V . Furthermore, limitations for the duty cycles d 1sc and d 2sc of the Split-Pi converter are considered:

0 ≤ d 1sc ≤ 1, 0 ≤ d 2sc ≤ 1.
Finally, an upper and lower limit for the UG is necessary as follows:

P min ug ≤ P ug (t) ≤ P max ug ,
with P min ug is -2100 W and P max ug is 4200 W . Since we aim the cost minimization, it is necessary to confine the power consumption from the UG, guaranteeing the full exploitation of the renewable sources and the energy storing devices.

Hierarchical constrained optimization-based control

In the previous sections, the system's characteristics, the control objectives and the different scenarios, which will be included in the optimization problem, were introduced. The combination of renewable sources, storing elements, variable consumers' demands and DC/DC converters lead to the multi-timescale dynamics. Hence, a hierarchical control problem divided in three levels will be proposed for the cost minimization with and without power losses in the central transmission network. At the high level, we deal with the slow dynamics, which, in this case, concludes to the power balancing in the central network. Then, at the lower level, the fast dynamics is introduced controlling the switches of the converters. The main goal is to reduce the electricity cost by minimizing the power generated by the UG, hence taking advantage of the PV power production and the ES system capacity. Considering all the above, the proposed hierarchical control approach is (Fig. middle level (battery scheduling): a tube-MPC tracking controller [START_REF] Langson | Robust model predictive control using tubes[END_REF]) is employed to mitigate the discrepancies among the reference and actual profiles under perturbation; low level (switching activity in the converter): an explicit control law for the duty cycles of the converter is provided to follow the tracking profiles obtained at the middle level.

In the following, each level is described in detail exploring two scenarios: i) the absence of power losses within the transmission lines; ii) the existence of power losses in the central transmission network.

High level control

The principal objective of the high level problem is the generation of reference trajectories in continuous time, exploiting the flat output representations provided in section 3. Therefore, an optimal scheduling for the battery charging and discharging is obtained and, in the meantime, the electrical power purchase from the UG is minimized. The general cost function considered at the high level is the following:

J = t f t0 ep(t)P ug dt, (4.3.1)
where ep(t) is the electricity price. Depending on the scenario, the form of the cost function changes, as it will be observed later in the detailed analysis.

Transmission line network without power losses

To determine the objective function at the high level for the reduced DC microgrid (Fig. 2.5.1) without power losses in the common DC bus, the dissipative elements will be erased from the power conservation equation (2.5.14). Therefore, considering that P R1 (t) = P R2 (t) = P R3 (t) = P R4 (t) = 0, (2.5.14) becomes:

P ug (t) + P pv (t)-P es (t) -P loads (t) = 0, (4.3.2a)
P ug (t) = P es (t) + P loads (t) -P pv (t) or P es (t) = P ug (t) -P loads (t) + P pv (t). (4.3.2b)

Hereinafter, we consider the minimization of the power coming from the UG while satisfying the ES system dynamics (2.5.1, 2.5.2a, 2.5.2b, 2.5.2c) and the constraints (4.2.17): The control variables of the optimization problem are the voltage, v b , and the current, i b of the battery. They are written in function of the flat output z(t) = q 2b (t) (3.3.26h) and, then, the B-spline parametrization is applied to obtain continuous-time optimal profiles for the charge, q 2b (t), hence, the voltage, v b , and the current, i b , of the battery. The main objective is the minimization of the electricity cost, ept, providing the optimal power generation from the ES, P es , and the UG, P ug , after considering the energy conservation equation (4.3.2a). Therefore, replacing in (4.3.3a) the power of the ES, P es , from (4.3.2b) with the control variables, a nonlinear optimization problem is obtained, since P es is in quadratic form.

min i b ,v b t f t0 ep(t)( P es (t) i b (t)v b (t) + P loads (t) -P pv (t))dt, (4.
In Fig. 4.3.1, the control scheme of the energy management problem, in the absence of power losses in the central transmission network, is presented. At each level, the control variables are illustrated. At the high level, the battery model and the constraints are presented in function of the B-splines. As aforementioned, the control variables are the voltage, v b , and the current, i b , of the battery. Optimal profiles are generated to use them at the middle level, where MPC tracking is applied and disturbances are considered. The control variable at this level is the output voltage of the Split-Pi converter, v sc out , which is linked to the voltage, v b , of the battery through Ohm's law. The discretized model of the battery is taken into account with its corresponding constraints in discrete-time. At this level, the main objective is the minimization of the discrepancies among the reference and the actual values. Afterwards, at the low level, the tracking profiles of the voltage, v b , and the current , i b , of the battery are applied to the ES PH system, as previously presented in section 2.5.1, in order to regulate the operation of the switches characterized by their duty cycles, d 1sc and d 2sc .

i b (t) = i R 1b (t), (4.3.4a) v b (t) = q 1b (t) C 1b , (4.3.4b)
i b (t) = i R 1b (t) = i R2b (t) + i C1b (t) = i 2b (t) + q1b = q1b (t) + q2b (t), (4.3.4c)
v b (t) = v 1b (t) = q 1b (t) C 1b , (4.3.4d)
where q 2b (t) considered to be the flat output of the battery as already introduced in sections 3.3.2.1 (3.3.26h) and also from the bicausal Bond graph in section 3.4.2.1 (3.4.17a). Therefore, considering also the flat representation of q 1b in (3.3.26g) and (3.4.17a):

q 1b (t) = R 2b C 1b ż(t) + C 1b C 2b z(t), (4.3.5a)
q 2b (t) = z(t), (4.3.5b)
we obtain the battery's voltage and current profile references expressed as:

v b (t) = 1 C 2b z(t) + R 2b ż(t), (4.3.6a) i b (t) = (1 + C 1b C 2b ) ż(t) + C 1b R 2b z(t). (4.3.6b)
Next in order is the B-spline approximation which ensures the continuous-time constraints validation. This will facilitate also the middle level controller to similarly respect them. Therefore, we replace z(t) with the B-splines approximation as in Definition 4.2.2:

z 3 (t) = N i=1 p i • B i,d (t) = PB d (t). ( 4
v b (t) = N i=1 1 C 2b P + R 2b P B (r) d (t), (4.3.8a) i b (t) = N i=1 1 + C 1b C 2b P + C 1b R 2b P B (r) d (t). (4.3.8b)
Then, employing the differentiation property of the B-splines (second property from section 4.2.2 (4.2.12)), the control points for voltage,v b , and the current, i b , are equal to:

p v b κi = 1 C 2b p i + R 2b (PM d ,d -1 S κ,d -1,d ) i , (4.3.9a)
p i b κi = 1 + C 1b C 2b (PM d ,d -1 S κ,d -1,d ) i + C 1b R 2b • (PM d ,d -2 S κ,d -2,d ) i . (4.3.9b)
Next, for v b (t) and i b (t) in (4.3.8a)-(4.3.8b), we have:

v b (t) = N i=1 1 C 2b p i + R 2b (PM d ,d -1 S κ,d -1,d ) i b i,d (t), (4.3.10a)
i b (t) = N i=1 1 + C 1b C 2b (PM d ,d -1 S κ,d -1,d ) i + C 1b R 2b • (PM d ,d -2 S κ,d -2,d ) i • b i,d (t), ∀t ∈ [τ κ , τ κ+1 ) . (4.3.10b)
Thus, the a priori optimization problem (4.3.3a-4.3.3f) is rewritten as follows:

min i b ,v b t f t0 P es (t) i b (t)v b (t)
+P loads (t) -P pv (t))dt =

v min,h b ≤ N i=1 p v b κi b i,d (t) ≤ v max,h b , (4.3.11c)
i min,h b ≤ N i=1 p i b κi b i,d (t) ≤ i max,h b , (4.3.11d)
q min,h 2b ≤ N i=1 p i b i,d (t) ≤ q max,h 2b , (4.3.11e)
P min,h ug ≤P ug ≤ P max,h ug . (4.3.11f)
The matrices S change across the knot sub-intervals [τ κ , τ κ+1 ]. Thus, constraints are also considered for each interval. For the P es constraint in (4.3.11f), we follow the power conservation equation (4.3.2a):

P min,h ug -P loads (t) + P pv (t) ≤ P es (t) ≤ P max,h ug -P loads (t) + P pv (t), (4.3.12)

where P es (t) = i b (t)v b (t) as aforementioned. Appendix B.1 demonstrates the detailed calculation of J es of (4.3.11a).

In this section, optimal profiles were obtained for the voltage, v b , and the current, i b , of the battery. Furthermore, the optimal power profiles for the ES and the UG were generated providing the optimal electricity cost. In the next section, we go further and consider in the optimization problem the power losses within the transmission network. The main idea is to find the optimal power routing while minimizing the power dissipation in the transmission lines and the electricity cost.

Transmission line network with power losses

In this section, the aforementioned flatness-based optimization approach is reformulated adding the power losses within the central network as in Fig. 4.3.4. Therefore, the dynamics of the system changes and the problem becomes more complicated because of the relations among the voltage, the current and the power of the DERs, the consumers' demand and the power losses. At this point, the dynamics of the central network, analyzed through its Bond graph in section 2.5.2.2 (see also Fig. 4.3.4), is considered. In Fig. 4.3.4, we observe that the notation on each connecting node has changed and, instead of v DC or i DC at the input of the ES, we have v es or i es . Since the power losses in the central transmission lines are considered, the voltage in the common DC-bus will not be stable. On each connecting node1 a different voltage value will appear, because of the voltage drop created by the resistors. Therefore, a different notation for the voltage and the current on each connecting node is necessary. In Fig. 4.3.3, the control scheme of the energy management problem, in the existence of power losses in the central transmission network, is illustrated. At the high level, the battery model and the constraints previously presented in section 4.3.1.1 in function of the B-splines are employed to the new objective function. The control variables are the voltage, v b , and the current, i b , of the battery, the voltages on the connecting nodes of the ES, v es , the PV, v pv , the UG, v u g, and the loads, v loads , and the α factor, which defines the relations among the duty cycles, d 1sc and d 2sc , in the Split-Pi converter, as already presented in section 2.4.2.1. The concept is similar to the previous case where the power losses in the common DC-bus are not considered. Optimal profiles are generated for the ES system, the duty cycles through the α factor and the voltages on the connecting nodes. The optimal power profiles for the ES, P es , and the UG, P ug , are provided minimizing the electricity cost and the power dissipation. At the high and the middle level, the whole dynamics of the ES system is taken into account, as it will be shown later, in order to link the power losses of the central transmission network to the battery dynamics. Afterwards, the optimal profiles obtained at the high level are used as references at the middle level for tracking under perturbation, where the control variable is the output voltage of the Split-Pi converter, v sc out . The discretized model of the battery is considered again and the main objective is the reduction of the deviations among the reference and the actual values. Finally, at the low level, the tracking profiles of the voltage and the current of the battery, as well as the input voltage of the ES system, v es , are applied to the ES PH model, as already described in section 2.5.1, to control the switching activity within the Split-Pi converter. At first, the power losses are taken into consideration in (2.5.14)

v es (t) = v b (t) + i b (t)R 1b α(t) + R 1sc α(t)i b (t) (4.3.20a) i es (t) = α(t)i b (t) (4.3.20b)
and the ES power is deduced as follows: Additionally , the following are also considered before writing the final form of the objective function of (4.3.15a): from (2.5.13), the voltage, v loads , on the connecting node n : 1 is written in function of the input voltage, v es , and the input current, i es , of the ES system:

P es (t) = v es (t)i es (t) = (v b (t) + i b (t)R 1b )i b (t) + R 1sc (α(t)i b (t))
v loads (t) = R 1 i es (t) + (1 + R 1 R 3 )v es (t) - R 1 R 3 v pv (t); (4.3.23)
from (2.3.4) and (2.5.17) the generated power from the PV system is written also in function of the voltages on the connecting nodes, v es , v pv , v ug : where N α is the number of control points p α i . Then, applying Theorem B.2.1 from Appendix B, the constraint of the α factor will be as follows:

P pv (t) = v pv (t)i pv (t) = v pv (t) v pv (t) -v es (t) R 3 - v ug (t) -v pv (t)
p α i > 1, ∀i = 1, . . . , N α . (4.3.28)
Since α lies in the interval of (0, 1), meaning that is positive, the following can be considered for the second part of (4.3.20a):

-

R 1sc |i b (t)| ≤ R 1sc α(t)i b (t) ≤ R 1sc |i b (t)|. ( 4
v b (t) + i b (t)R 1b α(t) -R 1sc |i b (t)| ≥ v min,h DC , (4.3.30a) v b (t) + i b (t)R 1b α(t) + R 1sc |i b (t)| ≤ v max,h DC . ( 4
p v b κ,i + R 1b p i b κ,i p α j -R 1sc |p i b κ,i | ≥ v min,h DC , (4.3.31a)
p v b κ,i + R 1b p i b κ,i p α j + R 1sc |p i b κ,i | ≤ v max,h DC , (4.3.31b)
where p v b κ,i and p i b κ,i are defined by (4.3.9) and κ, i and j ∈ N satisfy d -1 ≤ κ ≤ n -1 , κd + 2 ≤ i ≤ κ + 1 and 1 ≤ j ≤ N α . After calculation, which is given explicitly in Appendix B.3, the previous inequalities are rewritten as:

v min,h DC ≤ N i=1 Nα j=1 (p v b κ,i + R 1b p i b κ,i )p α j b i,j,d (t) - N i=1 Nα j=1 R 1sc |p i b κ,i |p α j b i,j,d (t) ≤ ≤ N i=1 Nα j=1 p v b κ,i + R 1b p i b κ,i -R 1sc |p i b κ,i | p α j b i,j,d (t), (4.3.32a) v max,h DC ≥ N i=1 Nα j=1 (p v b κ,i + R 1b p i b κ,i )p α j b i,j,d (t) + N i=1 Nα j=1 R 1sc |p i b κ,i |p α j b i,j,d (t) ≥ ≥ N i=1 Nα j=1 (p v b κ,i + R 1b p i b κ,i ) + R 1sc |p i b κ,i | p α j b i,j,d (t). (4.3.32b)
Concerning the constraint of the P ug (t), from (2.5.18), it is defined below:

P ug (t) = v ug (t)i ug (t) = v ug (t)[i R4 (t) -i R2 (t)] = v ug (t) v ug (t) -v pv (t) R 4 - v loads (t) -v ug (t) R 2 = = v 2 ug (t)( 1 R 4 + 1 R 2 ) -v ug (t) v pv (t) R 4 + v loads (t) R 2 . ( 4
P ug (t) = v 2 ug (t)( 1 R 4 + 1 R 2 ) - 1 R 4 v ug (t)v pv (t) -v ug (t) R 1 i es (t) + (1 + R 1 R 3 )v es (t) - R 1 R 3 v pv (t) . (4.3.34)
Including also (4.3.20a) and (4.3.20b), P ug (t) is denoted as:

P ug (t) = v 2 ug (t)( 1 R 4 + 1 R 2 ) + v ug (t)v pv (t)( R 1 R 2 R 3 - 1 R 4 ) - 1 R 2 v ug (t) R 1 α(t)i b (t) + (1 + R 1 R 3 ) v b (t) + R 1b i b (t) α(t) + R 1sc α(t)i b (t) = v 2 ug (t)( 1 R 4 + 1 R 2 ) + v ug (t)v pv (t)( R 1 R 2 R 3 - 1 R 4 ) - 1 R 2 v ug (t) (1 + R 1 R 3 ) v b (t) + R 1b i b (t) α(t) + (R 1 + R 1sc + R 1 R 1sc R 3 )α(t)i b (t) . (4.3.35)
Therefore, finally, the constraint (4.3.26b) is defined as:

P max,h ug ≥ v 2 ug (t)( 1 R 4 + 1 R 2 ) + v ug (t)v pv (t)( R 1 R 2 R 3 - 1 R 4 ) - 1 R 2 v ug (t) (1 + R 1 R 3 ) v b (t) + R 1b i b (t) α(t) -(R 1 + R 1sc + R 1 R 1sc R 3 )|i b (t)| , (4.3.36a)
P min,h ug ≤ v 2 ug (t)( 1 R 4 + 1 R 2 ) + v ug (t)v pv (t)( R 1 R 2 R 3 - 1 R 4 ) - 1 R 2 v ug (t) (1 + R 1 R 3 ) v b (t) + R 1b i b (t) α(t) + (R 1 + R 1sc + R 1 R 1sc R 3 )|i b (t)| . ( 4
(1 + R 1 R 3 ) v b (t) + R 1b i b (t) α(t) -(R 1 + R 1sc + R 1 R 1sc R 3 )|i b (t)| ≥ - R 2 v ug (t) P max,h ug -v 2 ug (t)( 1 R 4 + 1 R 2 ) -v ug (t)v pv (t)( R 1 R 2 R 3 - 1 R 4 ) , (4.3.37a)
(1 + R 1 R 3 ) v b (t) + R 1b i b (t) α(t) + (R 1 + R 1sc + R 1 R 1sc R 3 )|i b (t)| ≤ - R 2 v ug (t) P min,h ug -v 2 ug (t)( 1 R 4 + 1 R 2 ) -v ug (t)v pv (t)( R 1 R 2 R 3 - 1 R 4
) . Then, in the next subsection, the middle level of the hierarchical problem is described. Note that the reference profiles obtained at the high level for the battery current, i b , battery voltage, v b , and input voltage of the ES system v es will be denoted in the middle level as i ref b , v ref b and v ref es , respectively. Similarly for the α factor reference profile which is also mentioned at the middle level as α ref .

(1 + R 1 R 3 ) v b (t) + R 1b i b (t) α(t) ± (R 1 + R 1sc + R 1 R 1sc R 3 )|i b (t)| = (1 + R 1 R 3 ) N i=1 Nα j=1 (p v b κ,i + R 1b p i b κ,i )p α j b i,j,d (t) ± (R 1 + R 1sc + R 1 R 1sc R 3 ) N i=1 Nα j=1 p i b κ,i b i,j,d (t)
ep(t) Q cost (v b (t) + i b (t)R 1b )i b (t) + R 1sc (α(t)i b (t)) 2 -P pv (t) + P loads (t) + +Q loss ( v b (t) + i b (t)R 1b α(t) + R 1sc α(t)i b (t) -v loads (t)) 2 R 1 + (v loads (t) -v ug (t)) 2 R 2 + + (v pv (t) - v b (t) + i b (t)R 1b α(t) + R 1sc α(t)i b (t)) 2 R 3 + (v ug (t) -v pv (t)) 2 R
v min,h b ≤ N i=1 p v b κi b i,d (t) ≤ v max,h b , i min,h b ≤ N i=1 p i b κi b i,d (t) ≤ i max,h b , (4.3.39e) q min,h 2b ≤ N i=1 p i b i,d (t) ≤ q max,h 2b , ( 4

Middle level control

In here, a tube-MPC controller is introduced to track the output voltage reference profile, v ref sc out , of the Split-Pi converter under perturbation. Consider that the output voltage reference of the converter can be written in function of the battery current and voltage reference profiles obtained at the high level by solving (4.3.3a)-(4.3.3f):

v ref sc out (t) = v ref b (t) + i ref b (t)R 1b , (4.3.40)
according to Ohm's law obtained from the electrical circuit (Fig. 4.3.2). At this level, the system of the battery will be discretized through Euler explicit as presented below:

ẋ(t = kT s ) ≈ x(k + 1) -x(k) T s , (4.3.41)
where T s is the discretization sampling time. Consequently, using the Euler explicit method, the battery's dynamics is discretized having as state variables the charges of the battery, q 1b and q 2b , as input variable the output voltage from the Split-Pi converter, v sc out , and as output variables the current and voltage of the battery, v b and i b , denoted as:

x(k) = q1b (k) q2b (k) , ũ(k) = ṽsc out (k), ỹ(k) = ĩb (k) ṽb (k) , (4.3.42)
where we have the state vector x(k) ∈ R 2 , the input ũ(k) ∈ R and the output vector ỹ(k) ∈ R 2 . From the ES electrical circuit, we obtain some additional information according to which ĩb (k) = ĩsc (k) and ṽb (k) = q1b (k) C 1b .

Therefore, we obtain the discretized system as follows:

xb (k + 1) = Ax b (k) + B ũb (k), ỹ(k) = C xb (k) + Dũ b (k), (4.3.43) with A ∈ R 2×2 , B ∈ R 2 , C ∈ R 2×2 and D ∈ R 2 equal to: A =     1 - T s C 1b 1 R 1b + 1 R 2b T s C 2b R 2b T s C 1b R 2b 1 - T s C 2b R 2b     , B =    T s R 1b 0    , C =     - 1 C 1b R 1b 0 1 C 1b 0     , D =    1 R 1b 0    .

Transmission line network in the absence of power losses

Next a tracking MPC controller is proposed in which the cost penalizes the tracking error (the difference between the actual and the reference output profiles) over a finite prediction horizon N p : where Pes (k) = ĩb (k)ṽ b (k) as aforementioned. The above objective function is in quadratic form with nonlinear constraints and a variable electricity cost. Additionally, the profiles of the PV, the loads demand and the electricity price are employed as in Fig. 4.3.1.

min ũ(k) k+Np-1 i=k (ỹ(i) -ỹref (i)) Q ỹ (ỹ(i) -ỹref (i))+(ũ(i) -ũref (i)) R ũ ũ(i) -ũref (i) (4.

Error dynamics

At the high level, a desired profile is generated for the voltage v b and the current i b of the battery. The dynamics considered at the middle level has to follow these profiles (replaced by the output voltage of the Split-Pi converter v sc out ) in the best possible way. This approach is the so-called tube-M P C [START_REF] Langson | Robust model predictive control using tubes[END_REF] where an MPC law provides the nominal input (based on the nominal, noise-free dynamics) and the actual input adds to the nominal value a corrective term which counteracts the noise. Let us consider a nominal system ẋ(t) = f (x(t), ũ(t), 0) and an actual system ẋw (t) = f w (x w (t), ũw (t), w(t)). The error between them is denoted as s(t) = xw (t)x(t) and the error system is given as follows:

ṡ(t) = f w (x w (t), ũw (t), w(t))f (x(t), ũ(t), 0), (4.3.46) where x(t) ∈ X . Furthermore, a control signal is designed which is composed by a nominal input and a state feedback control as follows: ũw (t) = ũ(t) + K(x w (t), x(t)), (4.3.47) where ũ(t) ∈ Ũ and K(x w (t), x(t)) : X × X -→ R m , where min is the number of inputs. A set S ⊂ X ⊂ R n (n is the number of states) is called robust invariant for the error system in (4.3.46), if a feedback controller K exists for which ũ(t) + K(x w (t), x(t)) ⊆ Ũ ∈ R m such that for all = e(t 0 ) ∈ S and for all w ∈ W, s(t) ∈ S for t ≥ t 0 . S is the robust positively invariant (RPI) set of the error system (4.3.46) under a feedback control law K.

The tracking error, under certain assumptions, can be bounded by a RPI set. Since the profile to be tracked is generated at the high level, we can tighten the constraints considered in its design in order to guarantee reliability under noises (with the tightening factor being defined by the aforementioned RPI set). The pair of nominal input and nominal state (ũ b (k), xb (k) is generated by the application of the MPC law over the nominal dynamics of the battery (4.3.48b). Furthermore, the real dynamics (4.

Transmission line network with power losses

At this point, similarly to the case of the transmission central network without the power losses in section 4.3.2.1, we continue using the tube-MPC controller and the discrete dynamics of the system 4.3.43. However, since the power losses are taken into account, the MPC tracking problem is reformulated as follows: (4.3.52g)

min ũ(k) k+Np-1 i=k (ṽ es (i) -ṽref es (i)) Q ṽes (ṽ es (i) -ṽref es (i)) + (ũ(i) -ũref (i)) R ũ ũ(i) -ũref (i) (4.
In the last constraint P ug (t), the power losses must be considered in the power conservation equation as in (4.3.13a) in section 4.3.1.2. Therefore, (4.3.52g) is replaced by:

P min,m ug -Ploads (k) + Ppv (k) -PR1 (k) -PR2 (k) -PR3 (k) -PR4 (k) ≤ Pes (k), (4.3.53a) P max,m ug -Ploads (k) + Ppv (k) -PR1 (k) -PR2 (k) -PR3 (k) -PR4 (k) ≥ Pes (k), (4.3.53b)
where Pes (k) = ĩes (k)ṽ es (k). ĩes (k) and ṽes (k) are calculated with respect to the α factor reference profile, α ref , obtained at the high level as in (4.3.20a) and (4.3.20b):

v es (k) = v b (k) + i b (k)R 1b α ref,h (k) + R 1sc α ref (k), (4.3.54a) i es (k) = α ref (k)i b (k). (4.3.54b)

Low level control

In the low level control, the fast dynamics of the system is considered caused by the switching activity within the Split-Pi converter. Here, the duty cycles of the switches are supervised according to the tracking profiles we obtained in the middle level, that we will use as references.

In order to proceed to the analysis of the Split-Pi converter system, we follow the patent of United States Patent and Trademark Office No: US 6914420 B2 published on July 2005 [START_REF] Crocker | Power converter and method for power conversion[END_REF] as already described in section 2.2.1.3. The patent provides the relations among the input and output voltage, v sc in and v sc out , and the duty cycles, d 1sc and d 2sc , of the converter. As already mentioned, the Split-Pi always operates in down-conversion towards the positive direction and in up-conversion otherwise. This means that only the duty cycle d 2sc will function and change values from 0 to 1 (on/off switching between Sw 3 and Sw 4 as shown in section 2.2.1.3). Therefore, the duty cycle d 1sc is always equal to 0 (Sw 2 is always off and Sw 1 is always on as shown in section 2.2.1.3). Consequently, the control variable is d 2sc and through (2.4.20) and Ohm's law the following are deduced (also previously described in (4.3.18)

d 2sc (t) = 1 - v DC (t) -v 2 DC (t) -4(v sc out (t) -v b (t))(v sc out ) 2(v sc out (t) -v b (t)) . (4.3.59)
The last equation (4.3.59) can be valid only if v sc out (t) = v b (t).

In the case of the transmission network with power loss existence, the same procedure is followed at the low level since it concerns the internal supervision of the ES system.

Simulation results

In this section, the simulation results of the three levels of the hierarchical control design are presented. Table 4.4.1 depicts the parameters of the DC microgrid used for the simulations. To proceed, the reference profiles are taken into account, as presented in section 4.2.4.2. Furthermore, a collection of AGM 12-165 lead acid batteries is considered (165 Ah battery capacity) for the ES system. The DC microgrid is connected to the UG ( 4200 W maximum UG generation ) through a DC breaker, as shown also in Fig. 2.5.1. For the simulations we use MATLAB 2015a. Furthermore, YALMIP optimization toolbox [START_REF] Löfberg | Yalmip: A toolbox for modeling and optimization in matlab[END_REF] was chosen for both high and middle level which allows the use of the IPOPT solver [START_REF] Biegler | Large-scale nonlinear programming using ipopt: An integrating framework for enterprise-wide dynamic optimization[END_REF] capable to handle nonlinear optimization problems. For the low level, we design and implement the PH ES system (section 2.5.1) in MATLAB/Simulink in order to validate the proper operation of the switching activity within the Split-Pi converter.

Error dynamics

Firstly, the constraints of the current and the voltage of the battery should be adjusted at the high level according to the RPI set (4.3.50, 4.3.51) as presented before in section 4.3.2.2. At this point, we mention that disturbances are added at the nominal input variable of the discretized system (4.3.43, 4.3.2), the output voltage of the Split-Pi converter ṽsc out , equivalent to 5% of the difference between the minimum and the maximum value of the ṽref sc out in (4.3.44a). In Table 4.4.4, the constraints of the high level are depicted. These are confined according to the constraints chosen at the middle level and the defined RPI set below (4.4.1).

In

C 1b , C 2b 86400, 21600 [F ] R 1 , R 2 , R 3 , R 4 1 [Ω]
≤ v b (t) ≤ v max,h b 12 ≤ v b (t) ≤ 13 [V ] i min,h b ≤ i b ≤ i max,h b -9 ≤ i b (t) ≤ 9 [A] Constraints q min,m 1b ≤ q 1b (t) ≤ q max,m 1b
288.3 ≤ q 1b (t) ≤ 307.7

[Ah] q min,h 2b ≤ q 2b (t) ≤ q max,h 2b 72.5 ≤ q 2b (t) ≤ 77. PV, ES, in respect to the total generated power. The UG and and the PV source are producing the highest amount of the total power. In correlation with Fig. 4.4.3, in the case of the CU demand which is high during the day, the PV, generating 47% of the total energy, tries to satisfy the consumers demand and the rest is sold to the UG. While in the afternoon, both the UG and the ES contribute to the loads' supply especially after 4p.m. when the energy generated from the PV is decreasing. Overall, only 1% of the ES power is sold to the UG and 93% of the total energy produced is addressed to the consumers. On the other hand, for DU, when the electricity price is high during the day the demand is low. The spare of the PV power is sold to the UG, almost 12% of the total power consumed. Additionally the batteries are charged during the night, when the electricity costs less, and during the day, when the PV power is at its peak (9% for charging of the total power consumed). In general, we notice that in both cases whenever the PV power fully covers the consumer demand, the difference is stored to the ES system or sold to the UG. Additionally, the electricity cost was calculated for both cases being equal to 4.052 euros for the CU profile and 2.462 euros for the DU profile. The cost without using the ES system is equal to 4.737 euros for the CU profile and 2.644 euros for the DU profile, which clearly shows its importance. In general, the simulation results are as expected, where for the CU case, where the demand is high during the day, the PV mostly supports the consumers' demand. While, in the DU case, where the demand is high during the afternoon, the power generated by the PV is either stored to the ES or sold to the UG. The power stored is used later in the afternoon together with the UG to satisfy the consumers' demand. At this point, the lifetime of the battery is considered. For this reason, after carrying out several simulations for different number of control points, N , as in (4.3.11a)-(4.3.11f), we obtain Table 4.4.4 which shows the electricity cost, the computation time and the number of discharges of the battery for both load profiles. We observe that as the number of control points grows, the computation time and the number of battery discharges increase, whereas the electricity cost decreases. From the number of discharges, the lifetime of the battery can be defined. Although notice that a large number of discharges leads to a decrease of battery's capacitance. This results in a reduced life for the battery and, thus to higher operational costs (necessitated by its premature replacement). Therefore, the number of control points in the optimization problem needs to be kept as low as possible to extend the lifetime of the battery and increase the electricity cost savings. is, also, enough energy to sell to the UG (approximately 13% of the total power consumed). Hence, to make profit in a commercial environment, the use of larger installations for the renewable resources is important.

Time [h]

Likewise, the batteries' reference profiles for the current, voltage, state of charge are introduced where the constraints are verified. The electricity cost is also calculated which is 4.318 euros for the CUs and 2.713 euros for the DU. The cost without using the ES system is equal to 4.737 euros for the CU profile and 4.732 euros for the DU profile, where we observe that in the first case remains the same since the ES system usage is almost negligible. Although, in the case of the DU profile, the cost without the battery existence is a lot higher, since its usage is exploited in the best possible way. Afterwards, the power losses are depicted in Fig. 4.4.5 together with the constraint validation for the voltage on the four different nodes where the sources and the loads' systems stand. In the CU profile (Fig. 4.4.5a), the power losses are higher mostly because of the UG power purchase towards the loads during the day (transmission line R 2 : 0.24% of the total power loss), while in the DU is less (Fig. 4.4.5b). Furthermore, we observe that the power losses in lines R 3 and R 4 are similar for both cases because of the PV purchase towards the UG or the ES system. For transmission line R 1 , with the CU profile, the loss is around 0.12% of the total loss and higher than the DU case since a large amount of power is transmitted from the UG to the loads during the day. While, for the DU, the total loss of R 1 and R 2 , about 0.33%, mostly exists due to the fact that the load demand increases during the afternoon. It is visible from Table 4.4.5 that, finally, the power loss is higher in the CU scenario since also the total power produced is higher than in the case of the DU load profile. The calculation time of the simulation is 12 min for the DU load and 6 min for the CU load.

Time [h]

As previously, the simulations correspond to the expected results. Since, the CU is high during the day, considering also the power losses existing in the transmission line network, the UG is unable to charge the batteries, since it must generate power for the consumers in cooperation with the PV source. In the DU case, since the consumers' demand is low during the day, similar results are expected as before without including the power losses in the central transmission network. The PV power is able to charge the batteries and sell power to the UG. However, because of the power losses, the power weakens while passing through the transmission lines.

Simulation results of the middle level

This section introduces the results of the middle level using as reference the optimal profiles generated at the high level for the battery current, i b , and the voltage, v b as in Fig. , where the power produced and the power consumed are illustrated in respect to the total power, we observe that the optimal profiles obtained at the high level are very closely followed.

≤ ṽb (k) ≤ v max,m b 11.9 ≤ ṽb (k) ≤ 13.1 [V ] i min,m b ≤ ĩb (k) ≤ i max,m b -10.6 ≤ ĩb (k) ≤ 10.6 [A] Constraints q min,m 1b ≤ q1b (k) ≤ q max,m 1b 287.6 ≤ q1b (k) ≤ 308.4 [Ah] q min,m 2b ≤ q2b (k) ≤ q max,
In the transmission network without power losses (Table 4.4.7), for the CU profile, a slight difference of 1% in charging and discharging of the battery is observed as well as for the power production by the U G. The cost of the electricity increases at this case about 2% and from 4.052 raises to 4.140 for the real profile. The same situation is noticed also for the DU demand regarding the discharging of the battery, but since there is no difference in the UG power production, the electricity cost remains the same. Therefore, for both the CU and the DU profile, the reference profiles are closely followed. The calculation time of the simulation lasts around 3 min for both the CU and the DU demands. Similarly, for the transmission network with power losses, the same behavior is observed. There is small dissimilarities (less than 1%) between the UG power production/consumption reference and the real profiles as well as the battery's charging/discharging. The cost shortly changes, from 2.713 to 2.708 euros for the DU demand and from 4.318 to 4.308 euros . The electricity cost decreases slightly, approximately 0.2%, since the UG power production decreases and, in parallel, the power sold to the UG increases. Furthermore, the power loss is higher by 1% for the real profiles because of the slight raise of power generated and transmitted from and to the ES system. Although, in general the errors are low and the reference profiles are well followed. The simulations time endure around 14 min for the CU demand and 9 min for the DU demand.

Time [h]

Simulation results of the low level

In the following, we present the results obtained for the low level following the tracking profiles of the middle level for the battery current, i b , and the voltage, v b , under perturbation from Fig.

Comparisons with MPC for optimal profiles generation

As it has been previously described, the reference trajectories were obtained through differential flatness and B-spline parametrization. As a next step, we compare the reference trajectories generation obtained through differential flatness and B-spline parametrization with MPC, as presented also in [START_REF] Parisio | A model predictive control approach to microgrid operation optimization[END_REF] and [START_REF] Iovine | Power management for a dc microgrid integrating renewables and storages[END_REF]. For the simulation, the system with the power losses is considered (section 4.4.3). The simulation results are presented for the reference profiles obtained with MPC in Fig. 4.4.10 and Fig. 4.4.11. Similar trajectories are observed with a slight difference in the cost as it is depicted also in Table 4.4.9 of about 0.2 -3% for the DU and 2 -5% for the CU.

Because of the electricity cost difference between the optimal profiles obtained with MPC and flatness, we proceed to several simulations using MPC at the high level with different prediction horizons N p and sampling times T s . In Table 4.4.9, we notice that with MPC and a prediction horizon equal to 24 and a sampling time equal to 1200 s, more than 1000 s are needed to perform a simulation, approximately twice of the differential flatness calculation time for a 24h simulation. Furthermore, we observe that because of the discretization in MPC, the value of the sampling time greatly influences (for better or worse) the simulation's

Optimal profile generation of different scenarios

Scenario 1: (One week optimal profiles generation) In here, we perform a simulation for one week considering the system without power losses as in section 4.4.2. The simulation lasts approximately 3 h and is performed under N = 164 control points. In Fig. 4.4.12, the P ower Balancing and the constraints validation of the optimization problem are illustrated. For the reference profiles, the PV generated power holds the same profile every day, contrary to the load profile which varies and concerns a domestic environment. During the simulations, the battery is discharged 14 times (see also Fig. 4.4.12). Table 4.4.10 depicts the percentages of power produced and consumed by the sources and loads in respect to the total power. For the simulation, the same parameters as in Table 4.4.2 are considered. for the CU profile, the transmission lines R 1 and R 4 are with different values lower than the R 2 and R 3 considering a greater distance between the UG and the consumers, the PV and the ES than in the previous case; for the DU profile, the transmission lines R 1 and R 4 are equal to 1 Ω (as in section 4.4.3) and R 2 and R 3 are equal to 3 Ω keeping this way the renewable resource closer to the UG and the ES closer to the loads. In the first case, the electricity cost increases for about 30 cents compared to the results obtained in section 4.4.3 for the CU profile, since the power loss among the loads and the U G is higher and, as a result, the external grid generates more power to satisfy the consumers' demand. Furthermore, the PV gives priority to the battery. Since the power loss is less from the battery to the loads, the battery discharges and transmit power towards the consumers. Equal behavior is observed also for the second scenario, where similarly the external grid needs to generate more power because of the loss that exists in R 2 , while the profile of the ES power remains almost the same as in Fig. 4.4.4b. Therefore, from the aforementioned results, the number of PVs and batteries close to the load play an important role for the consumers' profit. The power generated from the renewable resources and the batteries must be able to overcome the power losses caused by the UG, reducing this way the UG power generation towards the microgrid. profiles are generated in case of faulted line events. We consider two scenarios with one faulted line each: R 2 = 0, the transmission line among the loads and the UG; R 4 = 0, the transmission line between the PV and the UG system.

Time [day]

0.5 [Ω] R 2 5 [Ω] R 3 5 [Ω] R 4 0.2 [Ω] Domestic R 1 , R 4 1 [Ω] R 2 , R 3 3 [Ω] Time [h]
In the first case, where the transmission line among the loads and the external grid doesn't work, the loss caused by the rest functional lines is higher than in the case of Fig. 4.4.4b. The UG generated power is distributed through three lines, hence the loss increases as does the electricity cost. Similar situation is observed in the case of line under fault between the PV and the UG system. The power loss is higher towards the UG system from the PV for selling power than in the case of Fig. 4.4.4b. This is because of the interruption of the direct connection among the UG and the PV and, as a result, the power passes through the remaining lines. In general, we observe in both cases, that, finally, the demand is fully satisfied. Therefore, the results validate the meshed topology of the network. In case of a line under fault, the optimization-based controller can meet the consumers' demand through the remaining transmission lines.

Conclusion

In this chapter, a multilevel supervision for a meshed DC microgrid was introduced. A constrained optimizationbased control approach was presented which solves the power balancing problem with or without power losses in the network. The three control levels were analyzed: the high, the middle and the low level. In a meshed topology, the optimization problem to solve becomes complicated since a model consists of multiple sources, different timescales, nonlinearities and constraint satisfaction at the same time. Therefore, we built a controller which could manage and take into consideration all the aforementioned factors. Afterwards, the reference profiles obtained by differential flatness and B-spline approximations were compared with the optimal profiles deduced by MPC proving that differential flatness represents an accurate and straightforward way to generate optimal profiles generation for power balancing optimization. Furthermore, considering different kind of scenarios, we analyzed the response of the system, verifying the operation of the meshed topology in the DC microgrid.

The main contributions of this chapter are summarized below:

the differential flatness and the associated B-splines approximations were used for optimal profile generation. The method handles successfully continuous-time constraints (validated not only at each sampling time, but for the entire sampling interval) and integral costs (along the entire simulation horizon); a multi-scale approach was implemented in the control architecture. Most of the works presented in the introduction of this chapter concentrate on a single aspect of microgrid control due to the large disparities in sampling times (from hours for load balancing to tenths of seconds for DC/DC converter switching). Here, each level is analyzed and an appropriate control law is given (which takes into account computation time limitations); overall, we may conclude that the PH implementation and the avoidance of discretization at the load balancing level provide significant improvements in comparison with standard approaches which discretize the dynamics and apply, e.g., MPC to solve the problem.

In the following chapter, the conclusions of this thesis are drawn.

Chapter 5

Conclusion and future research

Conclusions

In the present manuscript, a hierarchical control strategy under constraints was studied to generate the optimal power distribution within a meshed DC microgrid architecture (Fig. 2.2.1) composed by a renewable source (PV), an ES system, a collection of DC loads (e.g. EVs, LED lighting, printers, computers and the like) and DC/DC converters. To obtain its dynamical model, the DC microgrid was analyzed interpreting all its components as electrical circuits. Each component was studied separately, deriving its mathematical model from the associated Bond graph [START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF][START_REF] Karnopp | System dynamics: modeling, simulation, and control of mechatronic systems[END_REF] in chapter 2. From the Bond graphs, the PH state-space representation (van der Schaft et al., 2014) of every element was presented including the central transmission network. This approach was selected because it preserves the constitutive laws for the components together with the Kirchoff's law and the associated power continuity. Then, the hierarchical control approach was analyzed, divided into three levels, the high, the middle and the low level. An optimization-based approach under constraints was proposed at the high level under the combination of differential flatness and B-spline parametrizations for the PH model. Afterwards, at the middle level, a tracking MPC was developed to follow the optimal profiles obtained at the high level. Finally, the tracking profiles were used at the low level for the local supervision of the ES system. More precisely, at first, the PH state-space representation for the sources was introduced. The PH model of the solar cell and the KiBaM battery were derived from their Bond graphs. Afterwards, the DC/DC converters, the Split-Pi converter and the Ćuk converter, were analytically described and their PH representations were validated through simulations, considering as control variables their duty cycles. Next, the PH dynamical model of the overall central transmission network was provided replacing the transmission lines with RL circuits. The complexity of the resulted system led us to simplify the mathematical models according to the control objectives. In this work, the main goals were to develop an optimization problem capable to ensure power balancing in the central transmission network while simultaneously minimizing the power dissipation and the electricity cost. Therefore, high importance was given to the ES system and the central transmission-line network. A simplified version of the system was introduced, thus reducing the transmission lines model but still preserving the detailed dynamical model of the ES.

Before integrating the DC microgrid model into the optimization problem, a flatness-based representation [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF][START_REF] Levine | Analysis and control of nonlinear systems: A flatness-based approach[END_REF] was proposed which maintains the power-preserving interconnections and the power conservation properties of the PH system. Through the flatness representation, the states and inputs of the PH model were reformulated in function of the derived flat outputs. A significant advantage of this approach and its combination with the B-spline parametrization [START_REF] Suryawan | Constrained Trajectory Generation and Fault Tolerant Control Based on Differential Flatness and B-splines[END_REF] is the continuous-time constraints validation. Despite the fact that flat representation is an advantageous method for system analysis, the calculation of the flat outputs remains to date a challenging problem. In this work, the computation of the flat outputs was studied, primarily, implementing the algorithm for flat representation proposed in [START_REF] Franke | On the computation of flat outputs for nonlinear control systems[END_REF] for the PH models of the different components of the reduced microgrid model introduced in this work (Fig. 2.5.1). Through this algorithm the flat representation of each component was provided. In the sequel, the flat outputs resulting from the algorithm were compared with another set of flat outputs, proposed for each part of the system by selecting a number of states as flat outputs equal to the number of control inputs according to [START_REF] Levine | Analysis and control of nonlinear systems: A flatness-based approach[END_REF]. Comparisons among the flat output sets were realized to find the most appropriate set for the system's flat representation. From theory, it is known that the states and the control inputs are deduced through the flat outputs and a finite number of their derivatives. However, the simulations indicated that the more derivatives there are in the flat representation of the states 115 and inputs of the system, the more noise-affected the states and the control inputs are. Therefore, a suitable flat output set is the one requiring fewer derivatives.

Because of the complexity in the flat outputs' calculation, another method was investigated in order to develop a more straightforward way to derive the possible flat outputs for the PH systems. This method is called bicausality which was proposed to inverse the dynamics of the Bond graphs. Since flat representation and bicausality are both based on the notion of reversibility of the dynamics of the system, bicausality seems an appropriate tool to generate the possible flat outputs of the PH systems. The proposed method was elaborated through the electrical circuits of the battery and the Split-Pi converter, for which the flat outputs were found. Through simulations, the flat outputs were analyzed to decide which set is the best choice. the high level, where the flat representation of the ES system was integrated into the optimization problem. The flat model and its corresponding constraints were rewritten using B-spline parametrization for continuous-time constraint validation. Two cases were investigated: i) the cost minimization generating profiles for the utility grid (UG) and the energy storage (ES) power as well as the voltage and the current of the battery during charging and discharging; ii) the power dissipation was also added to the optimization problem cost to generate the optimal power distribution among the DC microgrid's components. Optimal profiles were produced for the UG and the ES power, the power losses, the voltages on the connecting nodes, the voltage and the current of the battery; the middle level, where a MPC tracking controller was developed to track the optimal profiles computed at the high level under perturbation. The error dynamics of the system was also analyzed. The objective function was introduced with and without power losses in the central transmission network; the low level, where the tracking profiles of the middle level were taken into account to be followed by the PH model of the ES system in order to supervise the switching activity inside the Split-Pi converter.

Furthermore, the high level flatness-based approach was compared with an economic MPC scheme. This comparison showed slight differences in the electricity cost among the power profiles generated by the flatnessbased approach and MPC. Noteworthy, the overall cost value was shown to depend on the sampling time and the length of the prediction horizon. This is the result of discretization which can affect the performance of the simulation. Other optimal profiles were generated for different scenarios such as in the case of fault events or for different distances among the components. In the case of a line under fault, it was proved that the meshed topology is reliable satisfying the consumers' demand even under faults. On the other hand, the relative distance between components was shown to be relevant. In particular, the renewable sources and the ES system benefit being close to the loads, thus avoiding power losses, ameliorating the operation of the microgrid, while also minimizing the electricity cost.

In the field of microgrids, many questions remain open. This work focused only on some particular modeling and control approaches. The control objectives include the power balancing of the DC microgrid and the scheduling of the energy storing elements. However, we firmly believe that the detailed PH models together with their flat representations and the hierarchical control approach developed in this work could further be used to cope with other challenging problems arising in this domain, as discussed in the following section.

Future improvements

Henceforward, the short term future work (Fig. 5.2.1) will be described. Since a modeling and an optimization approach were developed for the central transmission network, the operation of the system in case of unexpected events or power outages could be further studied. With the mathematical model developed for the central transmission network, fault mitigation is possible, improving the performance and the reliability of the system.

The reconfiguration of the system may be accomplished at the high level. The optimization problem developed including the power losses in the central transmission network could be used to predict the behavior of the system in case of unexpected events as has been already started in this work in section 4.4.7. Generating profiles in case of a line under fault, the corresponding line can be isolated by forecasting the new behavior of the remaining transmission lines. In such a way, the operation of the system continues taking into account the updated optimal profiles. Furthermore, these tools will help in storage sizing for functioning in islanded mode.

The stability of the system can be studied thorough Nonlinear Model Predictive Control (NMPC) which considers the PH models. Recall that, PH formalism is a useful approach to analyze the stability of the system, due to its properties related to power-preserving interconnections, the dissipative and the storing elements. According to (van der Schaft et al., 2014), the Hamiltonian, H, i.e. the energy stored in the physical systems, is regarded as a candidate Lyapunov function. Therefore, the stability of the PH system can be investigated over a MPC controller which handles constraints and nonlinearities in an improvement to other control methods used in PH systems, such as Energy-balancing passivity-based control or Interconnection and damping assignment passivity-based control [START_REF] Macchelli | Passivity-based control of implicit port-hamiltonian systems[END_REF]van der Schaft et al., 2014). For instance, a passivity-based controller could enforce locally the stability of the PH system and be used as a local controller for the terminal region used to formally prove the recursive feasibility in an NMPC framework [START_REF] Nguyen | On the use of a computed-torque control law for the terminal region of an nmpc scheme[END_REF]. Apart from this, at the middle level, the error dynamics can be further studied for the overall system, as was already started with the battery in section 4.3.2.2, defining the upper limits of the states.

In addition, the flat representation of the Split-Pi converter, presented in chapter 3, could be added, taking into account its input and output voltage with their flat representations and not only through the battery's flat model. This would be useful also to further analyze the switching activity within the converters. Optimal profiles for the states of the Split-Pi converter could be provided, considering power loss minimization. Further constraints could be considered for the states of the Split-Pi converter, thus ensuring a better regulation in the converters. Finally, the system could be developed in its primary form by adding back the transmission lines which were discarded in the reduced model. Not in the least, the transmission line model can consider inductors, relevant for the transitional behavior of the grid.

A variable sized storage component may also be considered. That is, we may consider electrical vehicles connected to the grid which temporarily can act as storage elements.

As a long term future work (Fig. 5.2.1), the bicausality notion, presented in section 3.4, should be further analyzed in order to develop a general methodology which could provide the flat outputs of the PH systems. In here, a first attempt was realized concerning the electrical circuits employed in this work and a first idea was developed on how we can proceed towards a more general solution which may cover not only the electrical circuits but also other physical systems. This information is essential and will facilitate the use of differential flatness in order to be able to exploit its advantages (e.g. preservation of the system properties) before incorporating the control part.

Another important issue which might be further analyzed is the correct topology and design of the proposed meshed DC microgrid to further reduce the power losses and the electricity cost, but also ameliorate the power quality. The number of batteries, the proper sizing, the number of PV panels, the type of renewable sources are some of the factors on which the performance of the system strongly depends. For example, the PV panels are a trustworthy choice to be part of a DC microgrid, since it is a DC source. However, when the solar irradiation is very low, its power generation is negligible. Therefore, other sources must be taken into consideration, e.g. air turbines, geothermal sources and the like.

A further aspect which may ameliorate the operation of the system is the distance among the sources and the consumers, which should be considered in the optimization problem, defined by the power loss in the transmission lines. Finally, the importance of the batteries' existence in such systems makes necessary the investigation of their lifetime. Improving the battery model, its lifetime could be considered, hence avoiding a fast deterioration of the energy storing element and, as a consequence, decreasing its cost and increasing the safe operation of the system. where Z [0] and z [1] are equal to: 7. Thus, the matrices P ⊥ 1,[0] and w [1] are presented below:

Z [0] =                -βq 2sc ηp 1sc 0 0 - βq 1sc q 2sc ηp 2 1sc 1 0 0 0                , ( A
P ⊥ 1,[0] =           0 0 1 0 0 0 0           , (A.

  4 (a) Output voltage of the Split-pi converter generated by the ES PH model. (b) Corresponding duty cycle d 2sc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.5.5 Total available charge of the KiBaM battery according to Fig. 2.5.4. . . . . . . . . . . . . . 37 2.5.6 Bond graph of the DC micogrid's central transmission network. . . . . . . . . . . . . . . . . 39 2.5.7 Bond graph of the DC micogrid's central transmission network with resistors for transmission lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.6.1 Electrical circuit of the reduced meshed DC micogrid. . . . . . . . . . . . . . . . . . . . . . 42 2.6.2 Central transmission network of the simplified meshed DC microgrid. . . . . . . . . . . . . . 43 3.1.1 General idea of chapter 3 where the differential flatness and the bicausality properties will be considered for a particular class of PH systems. . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.1 Electrical circuit of the ES presented in (3.3.18a)-(3.3.18g), where the input voltage v s is equivalent to the voltage entering from the central transmission network v DC . The resistor's current i R 1b equals to the battery's input current. . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3.2 Battery's input current i R 1b generated by the PH system in (3.3.19a)-(3.3.19g). . . . . . . . 53 3.3.3 Control inputs of the ES in function of the states generated considering the flat output set in (3.3.25a-3.3.25d). The signals in red depict the expected value of the control inputs. . . . 53 3.3.4 Control inputs of the ES in function of the second set of flat outputs as in (3.3.28). . . . . . 56 3.3.5 Electrical circuit of the Split-Pi converter presented in (3.3.31a)-(3.3.31e). . . . . . . . . . . 56 3.3.6 Output current of the Split-Pi converter generated by the PH model (3.3.31a-3.3.31e). . . . 58 3.3.7 Control inputs in function of the flat outputs generated by the algorithm for flat representation as in (3.3.38a)-(3.3.38d). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.3.8 Control inputs in function of the flat outputs as in (3.3.41). . . . . . . . . . . . . . . . . . . 60 3.4.1 Causality property of the Bond graph. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.4.2 Electrical RLC circuit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.4.3 Causal bond graph of the electrical RLC circuit depicted inFig.3.4.2. . . . . . . . . . . . . . 61 3.4.4 Bicausal assignment of the effort and flow sources. . . . . . . . . . . . . . . . . . . . . . . . 62 3.4.5 Relation among the ef pairs in the bicausal Bond graphs. . . . . . . . . . . . . . . . . . . 62 3.4.6 Definition of the shortest causal path between the system input and output. . . . . . . . . . 62 3.4.7 Assignment of the storing (C, I) and dissipative (R) elements. . . . . . . . . . . . . . . . . . 62 3.4.8 Final bicausal Bond graph of the electrical RLC circuit of Fig.3.4.2. . . . . . . . . . . . . . 63 3.4.9 Corresponding bicausal Bond graph with numbers on the arrows and the differential causalities. 63 3.4.10 Causal Bond graph of the KiBaM battery. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.4.11 Bicausal bond graph of the KiBaM battery. . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.4.12 Causal Bond graph of the Split-Pi converter. . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.4.13 Bicausal Bond graph of the Split-Pi converter. . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.4.14 Control inputs in function of the first flat output set in (3.4.28). . . . . . . . . . . . . . . . . 72 3.4.15 Control inputs in function of the second flat output set in (3.4.29). . . . . . . . . . . . . . . 72 4.2.1 B-spline curve (blue), its control polyline (dark grey) and convex hull (light grey) (Prodan et al., 2019). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.2.2 General MPC control scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.2.3 Profiles of electricity price, renewable power from the PV panel and loads for domestic and CU. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.3.1 Flowchart of the proposed hierarchical control strategy of the DC microgrid without considering the power losses in the central transmission network. . . . . . . . . . . . . . . . . . . . 84 4.3.2 Electrical network of the battery connected to the Split-Pi converter during charging mode, where v DC (t) = 400 V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.3.3 Flowchart of the proposed hierarchical control strategy of the DC microgrid considering the power losses in the central transmission network. . . . . . . . . . . . . . . . . . . . . . . . . 87 4.3.4 Central electrical network of the meshed DC microgrid system under the existence of power losses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.3.5 Electrical network of the battery connected to the Split-Pi converter during charging mode, where v es (t) varies according to the power losses in the central transmission network. . . . . 89 4.4.1 RPI sets for defining the constraints of high and middle level. . . . . . . . . . . . . . . . . . 97 4.4.2 Ultimate bounds for discrete time events of the model states, q1b and q2b following the CU load profile. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.3 (a) Power balancing, optimal reference profiles and state of charge of the ES system of the CU. (b) Power balancing, optimal reference profiles and state of charge of the ES system of the domestic load profile. The red lines at the upper and lower part represent the corresponding constraints. (Transmission line network without power losses) . . . . . . . . . . . . . . . . . 4.4.4 (a) Power balancing, optimal reference profiles and state of charge of the ES system of the commercial load profile. (b) Power balancing in the absence of power losses in the the transmission line network, optimal reference profiles and state of charge of the ES system of the domestic load profile. The red lines at the upper and lower part represent the corresponding constraints. (Transmission line network with power losses) . . . . . . . . . . . . . . . . . . . 4.4.5 (a) Power losses and voltage on the transmission network nodes for CU profile. (b) Power losses and voltage on the transmission network nodes for DU profile. The red lines at the upper and lower part represent the corresponding constraints. . . . . . . . . . . . . . . . . . 4.4.6 (a) Power balancing, tracking references, available charge and UG power of the commercial load profile. (b) Power balancing, tracking references, available charge and UG power of the domestic load profile. The red lines represent the corresponding constraints. (Transmissionline network without power losses) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.7 (a) Power balancing, tracking references, available charge and UG power of the commercial load profile. (b) Power balancing, tracking references, available charge and UG power of the domestic load profile. The red lines represent the corresponding constraints. (Transmissionline network with power losses) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.8 (a) Voltage and current tracking profiles for CU. (b) Voltage and current tracking profiles for DU. (Transmission-line network without power losses) . . . . . . . . . . . . . . . . . . . 4.4.9 (a) Voltage and current tracking profiles using the commercial load profile. (b) Voltage and current tracking profiles using the domestic load profile. (Transmission-line network with power losses) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.10 figure/Power balancing using commercial load profile. Comparisons of P es and P ug with optimal profiles obtained by MPC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.11 figure/Power balancing using domestic load profile. Comparisons of P es and P ug with optimal profiles obtained by MPC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.12 Power balancing and reference profiles within one week. . . . . . . . . . . . . . . . . . . . . 4.4.13 (a) Power balancing and power losses profiles for DU load demand. (b) Power balancing and power losses profiles for CU load demand. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.14 (a) Power balancing and power losses profile generation with R 2 equal to 0. (b) Power balancing and power losses profile generation with R 4 equal to 0. . . . . . . . . . . . . . . . 5.1.1 General hierarchical control scheme of the DC microgrid. . . . . . . . . . . . . . . . . . . . . 5.2.1 Short and long term future work. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1.1 Magnetic flux states of the inductors I 1sc and I 2sc of the Split-Pi converter in function of the flat outputs generated by the algorithm as in (3.3.25a)-(3.3.25d). The error between the magnetic flux states of (3.3.26a)-(3.3.26b) and (3.3.19a)-(3.3.19b) is depicted. . . . . . . . . A.1.2 Charge states of the capacitors C 1b and C 2b of the battery written in function of the flat outputs generated by the algorithm as in (3.3.25a)-(3.3.25d). The error between the charge states of (3.3.26g)-(3.3.26h) and (3.3.19f)-(3.3.19g) is depicted. . . . . . . . . . . . . . . . . A.1.3 Charge states of the capacitors C 1sc , C 2sc and C 3sc of the Split-Pi converter in function of the flat outputs generated by the algorithm as in (3.3.25a)-(3.3.25d). The error between the charge states of (3.3.26c), (3.3.26e), (3.3.26f) and (3.3.19c), (3.3.19d), (3.3.19e) is depicted. A.1.4 Magnetic flux states of the inductors I 1sc and I 2sc of the Split-Pi converter in function of the second set of flat outputs as in (3.3.28). The error between the magnetic flux states of (3.3.29a)-(3.3.29b) and (3.3.19a)-(3.3.19b) is depicted. . . . . . . . . . . . . . . . . . . . . . A.1.5 Charge states of the capacitors C 1sc , C 2sc and C 3sc of the Split-Pi converter in function of the flat outputs generated by the algorithm in (3.3.28). The error between the charge states of (3.3.29c), (3.3.29d), (3.3.29e) and of (3.3.19c), (3.3.19d) and (3.3.19e) is depicted. . . . . A.1.6 Charge states of the capacitors C 1b and C 2b of the battery written in function of the flat outputs generated by the algorithm in (3.3.28). The error between the charge states in function of the flat outputs (3.3.29f)-(3.3.29g)) and of (3.3.19f) and (3.3.19g) is depicted. . . A.1.7 (a) Signal response of the control inputs after filtering the derivatives shown in the ES system flat representation in (3.3.26i)-(3.3.26l). (b) Signal response of the control inputs after filtering the derivatives shown in the ES system flat representation with the set of random flat outputs (3.3.29h-3.3.29k) is depicted. . . . . . . . . . . . . . . . . . . . . . . . . 124 A.2.1 Magnetic flux states of the inductors I 1sc and I 2sc of the Split-Pi converter in function of the flat outputs generated by the algorithm as in (3.3.38a)-(3.3.38d). The error between the magnetic flux states of (3.3.39a)-(3.3.39b) and (3.3.32a)-(3.3.32b) is presented. . . . . . . . 127 A.2.2 Charge states of the capacitors C 1sc , C 2sc and C 3sc of the Split-Pi converter written in function of the flat outputs generated by the algorithm as in (3.3.38a)-(3.3.38d). The error between the charge states of (3.3.39c)-(3.3.39e)) and (3.3.32c)-(3.3.32e) is depicted. . . . . 127 A.2.3 Magnetic flux states of the inductors I 1sc and I 2sc of the Split-Pi converter in function of the second set of flat outputs as in (3.3.38). The error between the magnetic flux states of (3.3.42a)-(3.3.42b) and (3.3.32a)-(3.3.32b) is depicted. . . . . . . . . . . . . . . . . . . . . . 127 A.2.4 Charge states of the capacitors C 1sc , C 2sc and C 3sc of the Split-Pi converter in function of the flat outputs generated by the algorithm in (3.3.28). The error between the charge states of (3.3.29c), (3.3.42c), (3.3.42d), (3.3.42e) and of (3.3.32c), (3.3.32d), (3.3.32e) is depicted. . 128 A.2.5 (a) Signal response of the control inputs after filtering the derivatives shown in the Split-Pi converter system flat representation in (3.3.39f)-(3.3.39i). (b) Signal response of the control inputs after filtering the derivatives shown in the ES system flat representation with the set of random flat outputs (3.3.42f-3.3.42i) is depicted. . . . . . . . . . . . . . . . . . . . . . . . 128 Chapter 1

Figure 1 .0. 1 :

 11 Figure 1.0.1: General architecture of the microgrid.

Figure 1 .

 1 Figure 1.0.2: General structure of AC and DC microgrids.

Figure 1 . 1 . 1 :

 111 Figure 1.1.1: Hierarchical control general structure.

  Conferences and journals: Zafeiratou, I., Prodan, I. Boem, F., Lefèvre, L. Handling power losses in a DC microgrid through constrained optimization, In Proceedings of the 21st IFAC World Congress, Berlin, Germany, 2020 (submitted); Zafeiratou, I., Prodan, I., Lefèvre, L. On the flat representation for a class of port-Hamiltonian systems, In Proceedings of the 21st IFAC World Congress, Berlin, Germany, 2020 (submitted); Zafeiratou, I., Prodan, I., Lefèvre, L., Piétrac, L. (2019). Meshed DC microgrid hierachical control: a differential flatness approach, Journal of Electric Power Systems Research (accepted, December, 2019); Zafeiratou, I.., Prodan, I., Lefèvre, L., Piétrac, L. (2018). Dynamical modelling of a DC microgrid using a port-Hamiltonian formalism, In proceedings of the IFAC Conference on Mathematical Modeling, Vienna, Austria, 51(2), 469-474, DOI: 10.1016/j.ifacol.2018.03.079; Zafeiratou, I., Nguyen, D. V. A., Prodan, I., Lefèvre, L., Piétrac, L. (2018). Flatness-based hierarchical control of a meshed DC microgrid, In proceedings of the IFAC Conference on Nonlinear Model Predictive Control, Madison, Wisconsin (USA), 51(20), 222-227, DOI: 10.1016/j.ifacol.2018.11.017; Zafeiratou, I.., Prodan, I., Lefèvre, L., Piétrac, L. (2017). Control and optimization scheduling within a meshed DC microgrid, In IFAC 20th World Congress, Toulouse, France, HAL Id: hal-01656982. Presentations: Zafeiratou, I., Prodan, I., Lefèvre, L., Piétrac, L., Hierarchical control of a meshed DC microgrid under constraints, University College London, London, May 2019 (mobility project); Zafeiratou, I., Prodan, I., Lefèvre, L., Piétrac, L., Hierarchical control of a meshed DC microgrid under constraints, National Scientific Workshop -GT Micro-réseaux, Compiègne, France, November 2018; Zafeiratou, I., Nguyen, D. V. A., Prodan, I., Lefèvre, L., Piétrac, L., Hierarchical control of a meshed DC microgrid, National Scientific Workshop -GDR SEEDS et MACS, Smart Energies 2018, Paris, France, July 2018; Zafeiratou, I.., Prodan, I., Lefèvre, L., Piétrac, Modeling and load balancing of a DC microgrid using port-Hamiltonian formulation, National Scientific Workshop -Micro and Smart Grid, Paris, France, October 2017.

Figure 1 . 4 . 1 :

 141 Figure 1.4.1: Organization scheme of the thesis.

  Adibi et al. (2017), they use PH representation for voltage regulation in a small scale microgrid in islanded mode. Furthermore, Fiaz et al. (2013); Pham et al. (2015); van der Schaft et al. (2014) design and analyze a large scale power network composed by distributed generators (DGs) and loads. The authors use PH modeling to obtain a unified mathematical structure which contains the dynamics of the energy sources, the storing devices and their components and defines the energy preserving interconnections. Moreover, in electronics, they develop PH models for inverters, converters, transformers and other electrical devices for input/output voltage regulation, internal analysis of the components' dynamics and power balancing (Meshram et al. (2017), Bergna-Diaz et al. (2017), Schiffer et al. (

Figure 2 . 2 . 1 :

 221 Figure 2.2.1: Structure of the meshed DC microgrid.

Figure 2 . 2 . 2 :

 222 Figure 2.2.2: General circuit of a PV cell.

Figure 2 . 2 . 3 :

 223 Figure 2.2.3: Two-tank model of the KiBaM

Figure 2 . 2 . 4 :

 224 Figure 2.2.4: Electrical circuit model of the KiBaM in state of charge.

Figure 2 . 2 . 5 :

 225 Figure 2.2.5: Split-Pi converter electrical circuit.

Figure 2 . 2 . 6 :

 226 Figure 2.2.6: Two functionalities of the Spit-Pi converter. (a) and (b) represent the down-conversion, while (c) and (d) represent the up-conversion process in positive direction. Similarly, we have the same two functionalities at the opposite direction, where the load takes the place of the source and the source, v s , takes the place of the load.

Figure 2 . 2 . 7 :

 227 Figure 2.2.7: Electrical circuit of the Ćuk converter.

Figure 2 . 2 . 8 :

 228 Figure 2.2.8: Operation modes of the Ćuk converter.

2. 3

 3 Bond graph and port-Hamiltonian systems 2.3.1 Bond graph

  .3.1 (Karnopp et al., 2012; Paynter, 1961): energy-storing elements (inductors I, capacitors C), energy-dissipating elements (resistors R), energy sources (voltage S e and current S f sources).

Figure 2 . 3 . 1 :

 231 Figure 2.3.1: Bond graph (Junction structure).

Figure 2 . 3 . 2 :

 232 Figure 2.3.2: Electrical circuit.

Figure 2 . 3 . 3 :

 233 Figure 2.3.3: Bond graph of the electrical circuit of Fig. 2.3.2.

Figure 2 . 3 . 4 :For

 234 Figure 2.3.4: Causal strokes definition for effort and flow.

Figure 2 .

 2 Figure 2.3.5: (a) Integral causality of the storage elements. (b) Differential causality of the storage elements.

Figure 2 . 3 . 6 :

 236 Figure 2.3.6: Energy sources, storage and dissipative elements with their causal strokes structure.

Figure 2 . 3 . 7 :

 237 Figure 2.3.7: Causality assignment for junctions 0 and 1 in causal Bond graphs.

Figure 2 . 3 . 8 :

 238 Figure 2.3.8: Causality assignment on the sources.Figure 2.3.9: Causality assignment on the storing elements.

Figure 2 . 3 . 9 :

 239 Figure 2.3.8: Causality assignment on the sources.Figure 2.3.9: Causality assignment on the storing elements.

Figure 2 . 3 .

 23 Figure 2.3.10: Causality on the dissipative elements. Numbering of causal strokes.

Figure 2 . 3 .

 23 Figure 2.3.11: Routing elements: (a) Transformers. (b) Gyrators.

  2.2.1. Primarily, the Bond graph, as in Fig. 2.3.10, can be mathematically portrayed as a PH state-space representation. The explicit formulation of a PH system is written below (Duindam et al., 2009; van der Schaft et al., 2014):

Figure 2 .

 2 Figure 2.3.12: Dirac structure of a port-Hamiltonian system.

where e = e 1 e 2 e 3

 3 ... e n and f = f 1 f 2 f 3 ... f n . Definition 2.3.1 (Power-conservation property -van der Schaft et al. (2014)). : Consider a finite dimensional linear space F with E = F * . The subspace D, D ⊂ F × E is defined as a Dirac structure if and only if:

Definition 2 . 3 . 3 .

 233 (PH system -van der Schaft et al. (

Figure 2 . 4 . 1 :

 241 Figure 2.4.1: Bond graph of the meshed DC microgrid.

Figure 2 . 4 . 2 :

 242 Figure 2.4.2: Bond graph of the solar cell circuit.

Figure 2 . 4 . 3 :

 243 Figure 2.4.3: (a) I -V curve of the PV model. (b) P -V curve of the PV model.

Figure 2 . 4 . 4 :Figure 2 . 4 . 5 :

 244245 Figure 2.4.4: (a) Temperature profiles for three days of different months, on January, on March and on June, within a year. (b) Irradiation profiles corresponding to the days chosen for the temperature profiles in Fig. 2.4.4a.

Figure 2 . 4 . 6 :

 246 Figure 2.4.6: Bond graph representation of the KiBaM.

Figure 2 . 4 . 7 :

 247 Figure 2.4.7: Bond graph of the Split-Pi converter.

Figure 2 . 4 . 8 :

 248 Figure 2.4.8: Expected output voltage during downconversion in the Split-Pi converter.

Figure 2 .

 2 Figure 2.4.9: Corresponding duty cycle.

Figure 2 . 4 .

 24 Figure 2.4.10: Ćuk converter Bond graph.

Figure 2 . 4 .

 24 Figure 2.4.11: Expected output voltage.

Fig. 2 .

 2 4.11 illustrates the output voltage generated by the PH model of the Ćuk converter. The expected output voltage is well-followed. In Fig.2.4.12 the corresponding duty cycle, d cc , is depicted.

  .4.25) where v cc in = v s and v cc out = v R Lcc . The v s and v R Lcc variables can be found through the PH model in (2.4.21)-(2.4.22).

Figure 2 .

 2 Figure 2.4.13: Cyclic network of loads.

Figure 2 .

 2 Figure 2.5.1: Reduced architecture of the DC microgrid.

Figure 2 .

 2 Figure 2.5.2: Detailed electrical network of the KiBaM battery connected to the electrical circuit of the Split-Pi converter during charging mode.

Figure 2 . 5 . 3 :

 253 Figure 2.5.3: Corresponding Bond graph of the electrical network of the KiBaM battery connected to the Split-Pi converter.

Figure 2 . 5 . 4 :

 254 Figure 2.5.4: (a) Output voltage of the Split-pi converter generated by the ES PH model. (b) Corresponding duty cycle d 2sc . The battery's charging and discharging are controlled by the Split-Pi converter's switches. For the validation of the ES model, we consider several simulations. In Fig. 2.5.4 and Fig. 2.5.5 we provide the output

Figure 2 . 5 . 5 :

 255 Figure 2.5.5: Total available charge of the KiBaM battery according to Fig. 2.5.4.

.5. 4 )Fig. 2 .

 42 Fig. 2.5.6 addresses the corresponding Bond graph of Fig. 2.5.1, which provides also the PH model in (2.5.3). Since, the J c is equal to 0 (2.5.3) and the knowledge of the interconnections among the elements is indispensable, another method can be employed and is called kernel representation (van der Schaft et al., 2014): Definition 2.5.1. (Kernel representation -van der Schaft et al. (2014)): Every Dirac structure D ⊂ F × E, with E = F * can be described in kernel representation as follows: D = (f, e) ∈ F × E | Ee + Ff = 0, (2.5.5)

.5. 7 )

 7 Proposition 2.5.2. (Constrained input-output representation -van der Schaft et al. (2014)): Every Dirac structure, D ⊂ F × F * can be represented as

Proposition 2 . 5 . 3 .

 253 (Hybrid input-output representation -van der Schaft et al. (2014)): Consider that a Dirac structure is represented in kernel representation with matrices E and F , as mentioned in Theorem 2.5.1.

Figure 2 . 5 . 6 :

 256 Figure 2.5.6: Bond graph of the DC micogrid's central transmission network.

Figure 2 . 5 . 7 :

 257 Figure 2.5.7: Bond graph of the DC micogrid's central transmission network with resistors for transmission lines.

Figure 2 . 6 . 1 :

 261 Figure 2.6.1: Electrical circuit of the reduced meshed DC micogrid.

Figure 2 . 6 . 2 :

 262 Figure 2.6.2: Central transmission network of the simplified meshed DC microgrid.

  al., 2019; Nguyen et al., 2018; Poultney et al., 2018; Prodan et al., 2013) in order to validate the system's dynamics under constraints and generate optimal profiles for velocity, acceleration, forces or torques. Other works employ B-Splines parametrization of the flat outputs to ensure continuous-time constraints validation at the trajectory generation level Stoican et al. (2017). For instance, in Pham et al. (2015) the flat representation of an electro-mechanical elevator within a DC microgrid is taken into account to generate the elevator's optimal velocity while minimizing the energy dissipation. In Thounthong and Pierfederici (2010) the authors implement an algorithm which controls the operation of a DC power source (fuel cell generator) through a power converter. The flat representation of the converter is provided and a control law is developed to study the stability of the system. In[START_REF] Pahlevaninezhad | A new control approach based on the differential flatness theory for an ac/dc converter used in electric vehicles[END_REF] the flatness theory is considered in order to prove controllability of an AC/DC converter and write the overall dynamical model in function of the input power which consists also the control input and flat output of the system.

Figure 3 . 1 . 1 :

 311 Figure 3.1.1: General idea of chapter 3 where the differential flatness and the bicausality properties will be considered for a particular class of PH systems.

Definition 3 . 2 . 1 .

 321 (Differential flatness[START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF]): Consider the following nonlinear system:

n

  (t)), where N = 1, 2, .., n. Concerning the input matrix G ∈ R n×m in (3.2.3a), it is replaced by the matrix G M ∈ R n×m in (3.2.4) where n = m, i.e. the number of control inputs equals the number of flat outputs, hence the number of the corresponding states. The matrix G M ∈ R n×m , n = m, is constant, square and invertible. The same for matrix Q ∈ R n×n in (3.2.3a), which is replaced by the matrix Q M ∈ R n×n , n = m, since in (3.2.4) only the states representing the flat outputs of the system are included.

Theorem 3 . 3 . 1 .

 331 (Implicit function theorem (David C. Royster, 1998)): Consider a point (x 0 , y 0 )=(x 0 1 , x 0

  .3.10) At this point, the B [i] must have full column rank, otherwise another two matrices are added in (3.3.6):

Figure 3 . 3 . 1 :

 331 Figure 3.3.1: Electrical circuit of the ES presented in (3.3.18a)-(3.3.18g), where the input voltage v s is equivalent to the voltage entering from the central transmission network v DC . The resistor's current i R 1b equals to the battery's input current.

3. 3

 3 .2.1 ES flat output generation using the algorithm for flat representation from Section 3.3.1 Thus, as aforementioned in section 3.3.1, we start by eliminating the system's(3.3.19a-3.3.19g) inputs and we obtain its implicit representation as in(3.3.3), which contains nm equations written below:

  (3.3.20c) Afterwards, the linearization as in (3.3.4) of the implicit system is implemented. From equations (3.3.20a)-(3.3.20c), the following linearized system is derived: βq 2sc d q2sc + βdq 2sc q2sc + ηp 1sc d ṗ1sc + ηdp 1sc ṗ1scαηdp 1sc q 1scαηp 1sc dq 1sc + +µp 2sc d ṗ2sc + µdp 2sc ṗ2sc + γµdq 3sc p 2sc + γµq 3sc dp 2sc = 0, (3.3.21a) d q 1b + δψdq 1b + ψdq 2bµdp 2sc + d q3sc = 0, (3.3.21b) d q2bδψdq 1b + ψ dq 2b = 0. (3.3.21c)

Figure 3 . 3 . 2 :Figure 3 . 3 . 3 :

 332333 Figure 3.3.2: Battery's input current i R 1b generated by the PH system in (3.3.19a)-(3.3.19g).

  .3.3) less noise-affected data exist than in the control inputs generated by the secondly proposed set of flat outputs in (3.3.28) (Fig.3.3.4). This is because of the high number of derivatives that appears in the d 1sc , d 2sc and v s flat representation, in (3.3.30h), (3.3.30i) and (3.3.30j) respectively. For i R1b similar results are obtained since their values are equal for both flat representations (see (3.3.26l) and (3.3.29k)). As previously, in Appendix A.1.2 the simulations of the states in function of the second set of flat outputs (3.3.28) are also demonstrated. Additionally, the numerical derivatives are filtered to mitigate the noise and the simulation results are depicted in Fig.A.1.7b in Appendix A.1.2.

Figure 3 . 3 . 4 :

 334 Figure 3.3.4: Control inputs of the ES in function of the second set of flat outputs as in (3.3.28).

Figure 3 . 3 . 5 :

 335 Figure 3.3.5: Electrical circuit of the Split-Pi converter presented in (3.3.31a)-(3.3.31e).

Figure 3 . 3 . 6 :

 336 Figure 3.3.6: Output current of the Split-Pi converter generated by the PH model (3.3.31a-3.3.31e).

Figure 3 . 3 . 7 :

 337 Figure 3.3.7: Control inputs in function of the flat outputs generated by the algorithm for flat representation as in (3.3.38a)-(3.3.38d).

Figure 3 . 3 . 8 :

 338 Figure 3.3.8: Control inputs in function of the flat outputs as in (3.3.41).

Figure 3 . 4 . 1 :

 341 Figure 3.4.1: Causality property of the Bond graph.

Figure 3 . 4 . 2 :

 342 Figure 3.4.2: Electrical RLC circuit.

Figure 3 . 4 . 3 :

 343 Figure 3.4.3: Causal bond graph of the electrical RLC circuit depicted in Fig.3.4.2.

  .4.4, which is the opposite of the case of the causal Bond graph depicted in the previous chapter in Fig.2.3.8; finally, the differential causality is determined for the storing elements (C, I) and the arbitrary causality for the resistors (R) (Fig.3.4.7 and Fig.3.4.8 respectively). As previously mentioned in section 2.3.1, according to differential causality, the input's time derivative is equal to the ratio of the output as depicted in Fig.2.3.5b. Therefore, the capacitors obey to de dt the same relations as in the causal Bond graphs: e = f • R and f = e R .

Figure 3 . 4 . 4 :

 344 Figure 3.4.4: Bicausal assignment of the effort and flow sources.

Figure 3 . 4 . 5 :

 345 Figure 3.4.5: Relation among the ef pairs in the bicausal Bond graphs.

Figure 3 . 4 . 6 :

 346 Figure 3.4.6: Definition of the shortest causal path between the system input and output.

Figure 3 . 4 . 7 :Figure 3 . 4 . 8 :

 347348 Figure 3.4.7: Assignment of the storing (C, I) and dissipative (R) elements.

Figure 3 .

 3 Figure 3.4.9: Corresponding bicausal Bond graph with numbers on the arrows and the differential causalities. we conclude to the succeeding relations from the 0 and 1 junctions (Fig.3.4.9) 2 : f irst 1 junction f rom the lef t : e 1e 3e 4 = 0, (3.4.1a) f 2 = f 3 = f 4 ; (3.4.1b) 0 junction : e 4 = e 5 = e 6 , (3.4.1c) f 4f 5f 6 = 0; (3.4.1d) second 1phantom..junction : e 6 = e 7 , (3.4.1e)

From Fig. 3 .

 3 4.9, considering the aforementioned relations of the junctions (3.4.2) and (3.4.5), we can derive the states and inputs which represent the inverse dynamics of the RLC circuit (Fig.3.4.2):

. 4 .

 4 11 and Fig. 3.4.13) of the ES system (see the electrical circuit in Fig. 3.3.1), which includes the KiBaM battery and the Split-Pi converter, the relations among the efforts and flows are derived from the 0 and 1 junctions as in (3.4.1a)-(3.4.1f). Next, their kernel representations are formulated to prove that they compose a DS. Therefore, according to Definition 2.5.1, the relation Ee + F f = 0 must exist, which contains the storing elements, the dissipative elements and the external ports of the system: E e I1 e I2 . . . e en I e C1 e C2 . . . e Cn C e R1 e R2 . . . e Rn R u + (3.4.10)

  .4.11) D is the subspace D ⊂ F × E of the DS regarding Definition 2.3.2. Note that the capacitors are regarded as (f C , e C ) pairs, for which f C =ẋ and e C = ∇H(x), and the inductors as (f I , e I ) pairs, where e I =ẋ and f I = ∇H(x). The Hamiltonian, H, is considered in its linear form, i.e. the energy storage of the electrical circuit equals to H(x) = 1 2

Figure 3 . 4 .

 34 Figure 3.4.10: Causal Bond graph of the KiBaM battery.

  Figure 3.4.11: Bicausal bond graph of the KiBaM battery.

.4. 14 )

 14 Therefore, the bicausal Bond graph of the battery is also a PH system. By isolating the derivatives of the state on the left side of the equations in (3.4.13a)-(3.4.13c), we result in the primary PH state-space representation as in (3.4.12a)-(3.4.12a). The equation (3.4.14) gives one combination among the efforts and flows from the bicausal Bond graph. However, more combinations exist derived by the 0 and 1 junctions. Considering Definition 2.3.3, the capacitors are regarded as (f C , e C ), for which f C =ẋ and e C = ∇H(x), and the Hamiltonian, H, is considered in the linear case, i.e. the energy storage of the electrical circuit, equal to H(x) = 1 2

Figure 3 . 4 .

 34 Figure 3.4.12: Causal Bond graph of the Split-Pi converter.

Figure 3 . 4 .

 34 Figure 3.4.13: Bicausal Bond graph of the Split-Pi converter.

  the bicausal Bond graph of the Split-Pi converter is also a PH system. By placing the state derivatives on the left side of the equations in (3.4.20a)-(3.4.20e) (considering also the control inputs in (3.4.21a)-(3.4.21d)), we conclude to the primary PH state-space representation as in (3.4.19a)-(3.4.19e). Then, we consider likewise Definition 2.3.3, with the storing elements regarded as (f C , e C ) and (f I , e I ), where f C =ẋ and e C = ∇H(x) and e I =ẋ and f I = ∇H(x) respectively, and the Hamiltonian in the linear case (H(x) = 1 2

  .4.29) The flat representation of the Split-Pi converter is written below from (3.4.20a)-(3.4.20e), (3.4.21a)-(3.4.21d) and (3.4.24a) 3 :

Figure 3 . 4 . 15 :

 3415 Figure 3.4.15: Control inputs in function of the second flat output set in (3.4.29).

  al., 2015; Fallahzadeh-Abarghouei et al., 2018; Parisio et al., 2014; Vandoorn et al., 2013). Hereinafter, various methods used at each level of the hierarchical control architecture are introduced: the tertiary control (high level) deals with power flow and optimization and, in parallel, takes into consideration different economical aspects. At this level, the power flow is regulated among the different parts of the microgrid succeeding the load balancing within the transmission network and providing the most economical solution. Several methods for optimization have been developed including constraints, energy storage or uncertainties. Such a method was presented in Parisio et al. (2015); Prodan and Zio (2014), where they combine Model Preditive Control (MPC) with mixed integer programming for battery scheduling. In Farina et al. (2016); Velarde et al. (2017), the authors use MPC with chance constraints to deal with profile uncertainties. In Drgoňa et al. (

  . 4.2.1). The polynomials, connected together, constitute the control polyline of the B-spline curve as shown in Fig. 4.2.1.

Figure 4 . 2 . 1 :

 421 Figure 4.2.1: B-spline curve (blue), its control polyline (dark grey) and convex hull (light grey) (Prodan et al., 2019).

  b

Definition 4 . 2 . 2 .

 422 (B-spline parametrization of a flat output (Lyche et al., 2018)): The flat output z(t) is projected over N B-splines of order d as follows (Suryawan, 2012):

  al., 2017): P1. A B-spline curve of order d is C d -1 -continuous at any time instant included in the knot vector T and C ∞ -continuous at any other point (Stoican et al., 2017).

Figure 4 . 2 . 2 :

 422 Figure 4.2.2: General MPC control scheme.

Figure 4 . 2 . 3 :

 423 Figure 4.2.3: Profiles of electricity price, renewable power from the PV panel and loads for domestic and CU.

  4.3.1 and Fig. 4.3.3): high level (power flow optimization): optimal profiles for the battery current i b and voltage v b are provided to use at the middle level. Differential flatness and B-spline parametrization are combined for continuous-time constraint validation;

-Figure 4 . 3 . 1 :

 431 Figure 4.3.1: Flowchart of the proposed hierarchical control strategy of the DC microgrid without considering the power losses in the central transmission network.

Figure 4 . 3 . 2 :

 432 Figure 4.3.2: Electrical network of the battery connected to the Split-Pi converter during charging mode, where v DC (t) = 400 V .

  )(P loads (t) -P pv (t))dt(4.3.11a) subject to : the system dynamics (2.5.1) -(2.5.2b), (4.3.6a), (4.3.6b),(4.3.11b)

Figure 4 . 3 . 3 :

 433 Figure 4.3.3: Flowchart of the proposed hierarchical control strategy of the DC microgrid considering the power losses in the central transmission network.

Figure 4 . 3 . 4 :

 434 Figure 4.3.4: Central electrical network of the meshed DC microgrid system under the existence of power losses.

( 4 .

 4 3.37b) Next, through (B.3.1a) and (B.3.1b) in Appendix B.3, the left part of (4.3.37a) and (4.3.37b) is defined in function of the B-splines:

≤

  3.44a) subject to : the system dynamics (4.3.Pug (k) ≤ P max,m ug , (4.3.44f) with ỹref (k) = ĩref b (k) ṽref b (k) , the current and voltage references of the battery, and ũref (k) = ṽref sc out (k), the output voltage reference of the Split-Pi converter, taken with a sampling time, T s . The last constraint P ug (t), it is rewritten as shown below (4.3.2a): P min,m ug -Ploads (k) + Ppv (k) ≤ Pes (k) ≤ P max,m ug -Ploads (k) + Ppv (k), (4.3.45)

Definition 4 . 3 . 1 .

 431 (Robust positively invariant set (Yu et al., 2013)):

  3.48a) is affected by the bounded noise wb (k): xw b (k + 1) = Ax w b (k) + B ũw b (k) + wb (k), (4.3.48a) xb (k + 1) = Ax b (k) + B ũb (k), (4.3.48b) where wb (k) is the perturbation, xb (k) is the nominal state and xw b (k) is the real, noise-affected, state. Linking the nominal and actual inputs through the relation ũw b (k) = ũb (k) + K(x b (k)xw b (k)) allows us to write the tracking error dynamics: s b (k + 1) = (A + BK)s b (k) + wb (k), (4.3.49) where s b (k) = xw b (k)xb (k). For any pair (A, B) in (4.3.48b) controllable there exists a static feedback K such that (A+BK) is stable which means that there exists an RPI set S b for which s b (k) ∈ S b ∀ k ≥ k 0 holds. Such a set can be computed with the ultimate bounds method (Kofman et al., 2007) or iterative procedures (Olaru et al., 2010). Having s b (k) ∈ S b , it is equivalent with xw b (k) ∈ {x b (k) ⊕ S} (note that ⊕ is the Minkowski sum 2 ). In other words, the nominal xb (k) has to be chosen more conservatively than xw b (k). Thus, to ensure that xw b (k) ∈ Xb = {x min,m b ≤ xb (k) ≤ xmax,m b }, xb (k) must be confined as follows: xb (k) ∈ Xb S (4.3.50) Note that the restriction on x(k) translates to a similar restriction on ỹ(k): ỹb (k) ∈ Ỹb CS b , (4.3.51) where Ỹb is a shorthand notation for constraints (4.3.44c),(4.3.44d) and (4.3.44f). The symbol refers to the Pontryagin difference 3 . This tightening term CS b is considered in the profile generation at the high level at (4.3.10a)-(4.3.10b).

≤

  3.52a)subject to : the system dynamics (4.3.Pug (k) ≤ P max,m ug .

Fig. 4 .

 4 4.1 and Fig. 4.4.2 we present the RPI set S b and the nominal and noise-affected variables, respectively. At several time instants, the set was illustrated to highlight that the profiles are robust under

  (a) RPI set S b of the model states under perturbation. (b) RPI set of the battery's states constraints. (c) RPI set of the battery's voltage and current constraints.

Figure 4 . 4 . 1 :

 441 Figure 4.4.1: RPI sets for defining the constraints of high and middle level. bounded noise (i.e., the real trajectory lies in a tube centered around the nominal trajectory). The corresponding RPI set is given below, which actually depicts the set where the error among the nominal and actual values can be. Therefore, the error is bounded as follows (see also Fig. 4.4.1a): S b -0.38 [Ah] -0.11 [Ah] qw 1b -q1b qw 2b -q2b

Figure 4 . 4 . 2 :

 442 Figure 4.4.2: Ultimate bounds for discrete time events of the model states, q1b and q2b following the CU load profile. portrays the simulation for both the commercial use (CU) (Fig. 4.4.4a) and the domestic use (DU) (Fig. 4.4.4b) profiles. The figure showing the P ower Balancing (see Fig. 4.4.3a, 4.4.3b) contains the power generation coming from the ES and the U G within 24 hours taking into account the PV and the CU and DU consumers demand. For the CU, the demand is high during the day from 6a.m. to 4p.m.. On the other hand, the DU demand increases during the afternoon after 4p.m.. According to this two load profiles, the continuous-time reference trajectories for the P ug and the P es by solving the optimization problem expressed in function of the control points and B-splines without power losses (4.3.11a-4.3.11f). Afterwards, the reference profiles for the battery current, i b , voltage, v b , and charge, q 2b , are also generated in Fig. 4.4.3, where the constraints are satisfied. Note that the power positive sign indicates the power supplied to the microgrid.

Figure 4 . 4 . 3 :

 443 Figure 4.4.3: (a) Power balancing, optimal reference profiles and state of charge of the ES system of the CU. (b) Power balancing, optimal reference profiles and state of charge of the ES system of the domestic load profile. The red lines at the upper and lower part represent the corresponding constraints. (Transmission line network without power losses)

Figure 4 . 4 . 4 :

 444 Figure 4.4.4: (a) Power balancing, optimal reference profiles and state of charge of the ES system of the commercial load profile. (b) Power balancing in the absence of power losses in the the transmission line network, optimal reference profiles and state of charge of the ES system of the domestic load profile. The red lines at the upper and lower part represent the corresponding constraints. (Transmission line network with power losses)

Figure 4 . 4 . 5 :

 445 Figure 4.4.5: (a) Power losses and voltage on the transmission network nodes for CU profile. (b) Power losses and voltage on the transmission network nodes for DU profile. The red lines at the upper and lower part represent the corresponding constraints.

Figure 4 . 4 . 6 :

 446 Figure 4.4.6: (a) Power balancing, tracking references, available charge and UG power of the commercial load profile. (b) Power balancing, tracking references, available charge and UG power of the domestic load profile. The red lines represent the corresponding constraints. (Transmission-line network without power losses)

Figure 4 . 4 . 7 :

 447 Figure 4.4.7: (a) Power balancing, tracking references, available charge and UG power of the commercial load profile. (b) Power balancing, tracking references, available charge and UG power of the domestic load profile. The red lines represent the corresponding constraints. (Transmission-line network with power losses)

Figure 4 . 4 . 8 :

 448 Figure 4.4.8: (a) Voltage and current tracking profiles for CU. (b) Voltage and current tracking profiles for DU. (Transmission-line network without power losses)

Figure 4 . 4 .

 44 Figure 4.4.9: (a) Voltage and current tracking profiles using the commercial load profile. (b) Voltage and current tracking profiles using the domestic load profile. (Transmission-line network with power losses)

Figure 4 . 4 .

 44 Figure 4.4.12: Power balancing and reference profiles within one week.

Scenario 3 :

 3 (Faulted lines in the transmission network) In here, considering the DU profile, optimal

Figure 4 . 4 .

 44 Figure 4.4.13: (a) Power balancing and power losses profiles for DU load demand. (b) Power balancing and power losses profiles for CU load demand.

Figure 4 .

 4 Figure 4.4.14: (a) Power balancing and power losses profile generation with R 2 equal to 0. (b) Power balancing and power losses profile generation with R 4 equal to 0.

Figure 5 . 1 . 1 :

 511 Figure 5.1.1: General hierarchical control scheme of the DC microgrid.

Figure 5 . 2 . 1 :

 521 Figure 5.2.1: Short and long term future work.

4 . 6 .

 46 .1.9a) z[1] = dq 2sc .(A.1.9b)3. Afterwards, the matrices P + 1,[0] and P ⊥ 1,[0] as in (3.3.9) are described: Additionally, from (3.3.10), the matrices A [0] and B [0] are introduced below:A [0] = -aq 1sc p 1sc , (A.2.5a)B [0] = γµq 3sc + αµp 2sc q 1sc p 1sc αηp 1sc 0 -γµp 2sc . (A.2.5b)The matrix B [0] is different from 0 which proves the controllability of the system.5. Then, from(3.3.16) and (3.3.15), the u 1 and w 1 matrices are presented below:u [1]= ηp 1sc dp 1sc + µp 2sc dp 2sc + βq 2sc dq 2sc , Afterwards, since the rank of matrix B [0] is not of full column rank, we proceed as in Appendix A.1.1 with (3.3.11). Therefore, the matrices Z [0] and z[1] are equal to:

Figure A. 2 . 1 :Figure A. 2 . 2 :Figure A. 2 . 3 :

 212223 Figure A.2.1: Magnetic flux states of the inductors I 1sc and I 2sc of the Split-Pi converter in function of the flat outputs generated by the algorithm as in (3.3.38a)-(3.3.38d). The error between the magnetic flux states of (3.3.39a)-(3.3.39b) and (3.3.32a)-(3.3.32b) is presented.

Figure A. 2 . 4 :Figure A. 2 . 5 :

 2425 Figure A.2.4: Charge states of the capacitors C 1sc , C 2sc and C 3sc of the Split-Pi converter in function of the flat outputs generated by the algorithm in (3.3.28). The error between the charge states of (3.3.29c), (3.3.42c), (3.3.42d), (3.3.42e) and of (3.3.32c), (3.3.32d), (3.3.32e) is depicted.

  -r,d translation matrix from higher to lower degree basis functions M d ,d -r matrix performing the linear combinations of the lower-degree basis functionsT knot vector ∈ R N +d τ κ Sw 4 is 1d 1sc andSw 3 is d 1sc d 2sc duty cycle which for Sw 1 is 1d 2sc and Sw 2 is d 2sc Corresponding Bond graph of the electrical network of the KiBaM battery connected to the Split-Pi converter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.5.

	2.5.3
	Notation Description JS Junction structure S e voltage sources S f current sources F flow matrix E effort matrix f flow e effort SS source-sensor elements D e effort detector D f flow detector p magnetic flux q charge H the Hamiltonian, the energy stored in a physical system DS Dirac structure P H port-Hamiltonian x state vector of the system ∈ R Notation Description z flat output n number of states m number of control inputs l number of outputs N number of control points p n α number of control points p for the α factor d order of the B-spline d a order of the B-spline for the α factor r derivative of the B-spline κ number of knots t time in seconds Electrical network Notation Description AC alternative current DC direct current P CC point of common coupling P V, pv solar panel system ES, es energy storage system DOD depth of discharge U G, ug utility grid EV electrical vehicle system DER distributed energy resources DC/DC direct to direct current SN, sn smart node KiBaM Kinetic Battery Model P electrical power d duty cycle of the switches in the converters d 1sc duty cycle which for Sw switches of the Split-Pi converter e effort representing the voltage in an electrical circuit f flow representing the current in an electrical circuit q 1b available charge state of the KiBaM battery q 2b bound charge state of the KiBaM battery b KiBaM battery sc Split-Pi Converter cc Ćuk Converter c cyclic network s power source C capacitor I inductor R 1pv , R 2pv resistances of the PV cell Di diode R 1b resistance between the KiBaM battery and the Split-Pi converter R 2b resistance of the KiBaM battery R 1 , R 2 , R 3 , R 4 resistances of the transmission network R 1sc resistance among the DC network and the Split-Pi converters i current v voltage i sc out output current of the Split-Pi converter v sc out output voltage of the Split-Pi converter i sc in input current of the Split-Pi converter v sc in input voltage of the Split-Pi converter i DC current coming from the common DC-bus nc number of the connecting nodes in the electrical transmission network nl number of the transmission lines C n number of the capacitors p κ th knot I n number of the inductors

n u input vector of the system ∈ R m y output vector of the system ∈ R m J skew-symmetric matrix ∈ R n×n in function of time for linear case Q diagonal matrix of circuit's parameters ∈ R n×n R dissipation matrix of the system ∈ R n×n G control matrix of the system ∈ R n×m D throughput matrix of the system ∈ R m×m Differential flatness and B-splines parametrization i i th control point b i,d (t) i th B-spline of order d B(t) vector of the B-splines ∈ R d×N P vector of control points ∈ R 3×N vii S κ,d

Table 2 .

 2 2.1: Different phases of operation in Split-Pi converter.

	Modes of operation	Switching activity
		Positive direction
	Down-conversion	Phase 1:

Sw 2sc =Sw 4sc =1 and Sw 1sc =Sw 3sc =0 Phase 2: Sw 2sc =Sw 3sc =1 and Sw 1sc =Sw 4sc =0 Up-conversion Phase 1: Sw 2sc =Sw 4sc =1 and Sw 1sc =Sw 3sc =0 Phase 2: Sw 1sc =Sw 4sc =1 and Sw 2sc =Sw 3sc =0

  Sw 2sc =Sw 4sc =1 and Sw 1sc =Sw 3sc =0 Phase 2: Sw 1sc =Sw 4sc =1 and Sw 2sc =Sw 3sc =0 Sw 2sc =Sw 4sc =1 and Sw 1sc =Sw 3sc =0 Phase 2: Sw 2sc =Sw 3sc =1 and Sw 1sc =Sw 4sc =0

	Negative direction
	Down-conversion Phase 1: Up-conversion Phase 1: 2.2.1.4 Ćuk converter

Table 2 .

 2 3.1: Efforts and flows in several physical fields.

	Physical system	Effort e	Flow f
	Electrical	Voltage	Current
	Mechanical	Force Torque	Velocity Angular velocity
	Chemical	Chemical potential	Mass flow rate
	Thermal	Temperature Pressure	Entropy change rate Volume change rate

Table 2 .

 2 4.1: Parameters of the Split-Pi converter.

				Parameters	Values
				I 1sc , I 2sc	0.25 H
				C 1sc , C 2sc , C 3sc	0.0008 F
				R 1sc		0.1 Ω
	16	Output voltage of the Spplit-Pi converter
	14
	12			real value
				reference value
	10
	0	200	Time(s) 400 600	800	1000

Table 2 .

 2 4.2: Parameters of the Ćuk converter.

	Parameters	Values

  n is the number of resistors, I n is the number of inductors and nc is the number of connecting nodes. Therefore, according to (2.4.26), (2.4.27), (2.4.28), from the Bond graph in Fig.2.4.1, we conclude to the following PH state-space representation (where we have 7 connecting nodes and 7 RL circuits for the transmission lines):

	In	Rn	nc
	P In +	P Rn +	P nc = 0,	(2.4.28)
	i=1	i=1	i=1
	where R

.

4.27)

where nl is the number of the transmission lines. Additionally to the PH model, we always consider the conservation property inside the transmission network, as written below:

Table 2

 2 .4.1 plus the values for C 1b and C 2b equal to 86400F and 21600F respectively. The capacitors values have been calculated are based on the AGM 12-165 battery parameters (Victron Energy, 2015).

  ). Proposition 3.2.2. The particular class of PH systems considered in this work (sections 2.4.1.1, 2.4.1.2, 2.4.2, 2.4.3, 2.5.1, 2.5.2.2) are differentially flat.

  .3.12) If the B [i] matrix is equal to 0 in (3.3.10), w [i+1] coordinates disappear and (3.3.10) is replaced by a differential equation which corresponds to a non controllable system.

	5. As a next step, the variables w [i+1] are eliminated from (3.3.10), using the orthogonal complement of
	B [i] for which B ⊥ [i] B [i] = 0. Therefore, (3.3.10) becomes:

  .3.14) 6. In addition, w [i+1] is determined for every step of the algorithm, multiplying in (3.3.10) by the pseudoinverse of the B [i] , B + [i]

  .3.15) Concerning u [i+1] , multiplying (3.3.7) with P 1,[i] , every next step of u [i+1] coordinates can be computed in function of u [i] as follows:

  for the renewable source, 28 PV modules of DS-100M PV model (180 W peak PV generation of 72 cells) are taken into account. The PV power profile were obtained after simulations based on the PV model presented in section 2.4.1.1. External temperature and irradiation data were gathered for a whole day in June (CIAT company, 2014). A profile with high power generation is selected to observe the controller response regarding the power sold to the UG or stored to the ES system; for the consumer's demand, real historical data gathered from National Renewable Energy Laboratory (2016) and Department of Energy Office of Energy Efficiency & Renewable Energy (2017) are used, one for commercial usage (CU) (4308 W peak demand) (National Renewable Energy Laboratory, 2016), where the demand is higher during the day, and one for domestic usage (DU) (3901 W peak demand) (Department of Energy Office of Energy Efficiency & Renewable Energy, 2017), where the demand increases after 4p.m. in the afternoon.

	Price [euros/kWh]	0.12 0.14	0	2	4	6 Electricity price reference profile 8 10 12 14 16 18	20	22	24
										Time [h]
		6000				PV panel power reference profile
	[W]	4000
	pv	2000
	P		0
				0	2	4	6	8	10	12	14	16	18	20	22	24
										Time [h]
						Consumers demand reference profile for domestic use
	[W]	4000
	P load	2000
				0	2	4	6	8	10	12	14	16	18	20	22	24
										Time [h]
	[W]	4000
	P load	2000
				0	2	4	6	8	10	12	14	16	18	20	22	24
										Time [h]

variable cost is considered: 0.147 [euros/kW h] from 4p.m. to 10p.m. and 0.116 [euros/kW h] for the rest of the day, with equal selling and buying values for the UG;

Table 4 .

 4 2.1: Number of discharge cycles of the battery.

	Number of cycles	DOD
	400	20%
	600	50%
	1500	70%

Table 4 .

 4 2.2: Boundaries of the battery's variables.

  .3.7) with control points P = [p 1 ...p N ] and basic functions of order d B d (t) = [B 1,d (t)...B N,d (t)]. From (4.3.6a) and (4.3.6b), the control variables are expressed in function of the B-splines as follows:

  and the P ug (t) is equal to: P ug (t) = P es (t) + P loads (t)-P pv (t) + P es (t) = v es (t)i es (t) is the input power on the third connecting node, P ug (t) = v ug (t)i ug (t) is the output power on the second connecting node, P pv (t) = v pv (t)i pv (t) is the output power on the fourth connecting node and P loads (t) = v loads (t)i loads (t) is the input power on the first connecting node of the central transmission network (see also Fig.4.3.4 and Fig.4.3.2). P R (t) = P R1 (t) + P R2 (t) + P R3 (t) + P R4 (t) is the total power loss created by the four resistors existing in the central network.Hereinafter,the new objective function is analytically calculated, which will be in the following form: system and the input voltage, v b (t), and current, i b (t), of the battery are obtained below:v sc in (t) = v es (t) -R 1sc i es (t)(4.3.19a)

	of the ES v sc out (t) = v b (t) + i b (t)R 1b	(4.3.19b)
						i sc in (t) = i es (t)	(4.3.19c)
					i sc out (t) = i b (t)	(4.3.19d)
	Therefore, combining the above equations in (4.3.19) and in (4.3.18), we conclude to the relations below:
				t f
			J =	t0	ep(t) P es (t) + P loads (t) -P pv (t) + P R (t) dt	(4.3.14)
	and the restored optimization problem will be:
		t f
	min i b ,v b ,α,ves,vpv,vug	t0	ep(t) P es (t)	(t)+P loads (t) -P pv (t) + P R (t) dt,	(4.3.15a)
			i b (t)v b (t)
		subject to : the system dynamics (section2.5.2.2)and(2.5.1) -(2.5.2b), the power conservation (4.3.13a),	(4.3.15b) (4.3.15c)

R (t) or (4.3.13a) P es (t) = P ug (t) -P loads (t) + P pv (t) -P R (t),

(4.3.13b)

where P

  es (t), i es (t), P es (t) are defined by (4.3.20a), (4.3.20b), (4.3.21) respectively. Since the variables v es (t) and i es (t) are written in function of the voltage and the current of the battery, v b and i b , they can be written also in function of the B-splines according to (4.3.10a) and (4.3.10b).

								2	(4.3.21)
	Substituting (4.3.16) in the cost function in (4.3.14), the cost function becomes:
	J =	t f	ep(t) Q cost P es (t)	-P pv (t) + P loads (t) +Q loss	(v es (t) -v loads (t)) 2 R 1	+	(v loads (t) -v ug (t)) 2 R 2
		t0	ves(t)ies(t)
				+	(v ug (t) -v pv (t)) 2 R 4	+	(v pv (t) -v es (t)) 2 R 3	dt.
								(4.3.22a)
	where v

  where loads is the soft constraint to relax the load's demand in order to ensure the feasibility of the optimization problem.The voltages on the connecting nodes, v es , v ug , v loads and v pv , obey to (4.2.18). The v loads and i loads are also written in function of the B-splines, since they are in function of v es and i es as in (4.3.23).

	To ensure constraint validation in continuous-time, (4.3.26a) and (4.3.26b) must be deduced in function of the B-splines. The α factor is, firstly, parametrized considering 1 α(t) ∈ (1, +∞) and a set of B-splines basis
	functions of order d α as follows:
			1 α(t)	=	Nα j=1	p α j b i,d α (t)	(4.3.27)
								R 4	;	(4.3.24)
	from (2.3.4) and (2.5.16), the constraint for the consumer's demand is considered as follows:
	P loads (t) -loads ≤ v loads (t)	v es (t) -v loads (t) R 1	-	v loads (t) -v ug (t) R 2	≤ P loads (t) + loads ,	(4.3.25)
								Two last
	constraints remain to be found for the α(t) factor and the P ug (t) in (4.3.15h) (taking into account also
	(4.3.13a)):
				α(t) -1 > 1,	(4.3.26a)
	P min,h	P min,h ug	ug ≤ P ug (t) ≤ P max,h	.	ug	.	(4.3.26b) (4.3.26c)

ug

≤ P es (t) -P pv (t) + P loads (t) + P R1 (t) + P R2 (t) + P R3 (t) + P R4 (t) ≤ P max,h

Table 4 .

 4 4.1: Parameters of the variables.

	Variable	Values	Units
	R 1sc , R 1b , R 2b	1, 0.025, 0.088	[Ω]
	I 1sc , I 2sc	0.25, 0.25	[H]
	C 1sc , C 2sc , C 3sc	0.0008, 0.0008, 0.0008	[F ]

Table 4 .

 4 4.2: Variables and constraints for the high level.

		Variable	Values	Units
		Without power losses
		N as in (4.3.11a)	18
	High level	d as in (4.3.11a)	4
		With power losses
		N as in (4.3.32a) (4.3.32b), (4.3.38)	27
		N a as in (4.3.32a) (4.3.32b), (4.3.38)	18
		d=d a as in (4.3.11a)	4
		v min,h b

Table 4 .

 4 4.3 presents the percentages of the power produced or consumed by the sources and loads, UG,

	5	[Ah]

Table 4 .

 4 4.3: Percentage of power in respect to the total power produced or consumed without power losses in the central network.

	Load profile	Power	Power produced [%] Power consumed [%] Electricity cost [euros]
		P ug	46%	1% sold to the UG	4.052
	Commercial	P es P pv	7% 47%	6% for ES charging -	--
		P loads	-	93% for load usage	-
		P ug	40%	12% sold to the UG	2.462
	Domestic	P es P pv	9% 51%	9% for ES charging -	--
		P loads	-	79% for load usage	-

Table 4 .

 4 4.4: Results for different numbers of control points.

	Load profile	N as in (4.3.11a)	Electricity cost [euros]	Computation time [s]	ES discharges
		18	4.052	157	2
		27	3.991	268	3
	Commercial	36	3.576	525	7
		45	3.409	865	9
		54	3.188	1333	11
		18	2.462	136	2
		27	2.505	248	2
	Domestic	36	2.193	482	7
		45	2.002	699	8
		54	1.797	916	9

Table 4 .

 4 4.4.3. Firstly, in Table4.4.6, the parameters are set for the middle level simulations. As previously mentioned, in the middle level, we use MPC for reference tracking with a prediction horizon N p equal to 5 and a sampling time T s equal to 5 Additionally, a perturbation is added to the system that is always lower than the difference between the maximum and minimum state value defined previously by the RPI set.In Fig.4.4.6 and 4.4.7, we observe the tracking profiles of the P ower Balancing, and the control input, v sc out (the output voltage of the Split-Pi converter), which is in function of the current, i b , and the voltage, 4.6: Variables and constraints for the middle level.

	min. The system, at this point, is discretized as in (4.3.48a) and (4.3.48b) and it follows the optimization
	problem under constraints as in (4.3.44a)-(4.3.44f) without losses or as in (4.3.52a)-(4.3.52g) with losses in
	the transmission network.

  b of the battery as in(4.3.40). From the figures and the tables, Table4.4.7 and Table 4.4.8

	additional constraint (with power losses)	m ≤ Pug (k) ≤ P max,m 2b ug ≤ v ug,pv,es,loads (t) ≤ v max,h P min,m ug v min,h DC DC	72.3 ≤ q2b (k) ≤ 77.7 -2100 ≤ Pug (k) ≤ 4200 370 ≤ v ug,pv,es,loads (t) ≤ 430	[Ah] [W ] [V ]

v

Table 4 .

 4 4.7: Percentage of power in respect to the total power produced by the sources or consumed by the sources and the loads without power losses in the transmission lines.

					Power	Power
	Load profile Power	Power produced [%]	Power consumed [%]	production difference from	consumption difference with
					high level [%]	high level [%]
		P ug	47%	1% sold to the UG	1%	0%
	Commercial	P es P pv	6% 47%	7% for ES charging -	-1% 0%	1% -
		P loads	-	92% for load usage	-	-1%
		P ug	40%	12% sold to the UG	0%	0%
	Domestic	P es	9%	11% for ES charging	0%	-2%
		P pv	51%	-	0%	-
		P loads	-	77% for load usage	-	-2%

Table 4 .

 4 4.8: Percentage of power in respect to the total power produced or consumed with power losses in the central transmission network.

						Power	Power
	Load profile Power	Power produced [%]	Power consumed [%]	production difference from	consumption difference with
						high level [%]	high level [%]
			P ug	49.24%	1.85% sold to the UG	-0.62%	0.07%
	Commercial	P es	2.08%	1.84% for ES charging	1.29%	1.05%
			P pv	48.68%	-	-0.62%	-
			P loads	-	94.74% for load usage	-	0.75%
			P loss	-	1.57%	-	1.01%
			P ug	41.27%	13.15% sold to the UG	-0.87%	0.08%
	Domestic	P es	7.95%	7.73% for ES charging	1.37%	0.97%
			P pv	50.78%	-	-0.5%	-
			P loads	-	77.71% for load usage	-	-1.05%
			P loss	-	1.41%	-	0.77%
			Output voltage of the battery
	b [V] 12.6
	Voltage V	12 12.2 12.4 11.8		reference profile real profile
		0 2 4 6 8 10 12 14 16 18 20 22 24
				Time [h]
				Time [h]

  4.3.52a and Fig. 4.4.7. For both central transmission networks with or without power losses, the ES PH model (section 2.5.1) is considered, designed in MATLAB/Simulink. Afterwards, in the case of the powerlossless transmission network, the current, i b , and the voltage, v b , of the battery are considered as the reference profiles to follow, taking into account the control law developed in section 4.3.3 (4.3.59). The voltage coming from the central transmission network is stable equal to 400 V . For the transmission network with power losses, apart from the voltage and current reference profiles of the battery, the reference profile of the voltage, v es , entering from the central transmission network is also considered.

				Output voltage of the battery	Output voltage of the battery
	Voltage V b [V] 12.2 12 12.1	Time [h] 0 2 4 6 8 10 12 14 16 18 20 22 24 reference profile real profile	b [V] 12.6 12 12.2 12.4 Voltage V 11.8	Time [h] 0 2 4 6 8 10 12 14 16 18 20 22 24 reference profile real profile
		4		Output current of the battery
	Current I b [A]	-2 0 2		reference profile real profile	Current I b [A]
		-4	0 2 4 6 8 10 12 14 16 18 20 22 24	0 2 4 6 8 10 12 14 16 18 20 22 24
				Time [h]	Time[h]
				(a)

Table 4 .

 4 4.9: Simulation results obtained for optimal profiles with MPC.

	Load	Prediction horizon N p [h]	Sampling time T s [s]	Electricity cost [euros]	Power loss [%]	Calculation time [s]
		24	1800	4.569	1.72%	1094
		24	1200	4.433	1.73%	522
		15	1800	4.418	1.74%	262
	Commercial	15	1200	4.377	1.74%	251
		10	1800	4.441	1.74%	329
		10	1200	4.361	1.76%	310
		10	600	4.379	1.77%	763
		24	1200	2.801	1.63%	1315
		15	1800	2.747	1.65%	475
	Domestic	15 15	1200 600	2.763 2.719	1.64% 1.63%	400 751
		10	1800	2.788	1.64%	420
		10	1200	2.734	1.65%	452
	off-line simulations can predict system's behavior for different scenarios as in the case of lines under
	fault or long-time period simulations.

Table 4 .

 4 4.10: Simulation results of Scenario 1.

	Power	Power produced [%]	Power consumed [%]	Electricity cost [euros]
	P ug	43.5%	13.4% sold to the UG	21.169
	P es	8.2%	8.2% for ES charging	-
	P loads	-	78.4% for load usage	-
	P pv	48.3%	-	-
	Scenario 2: (Different values on the central lines resistors) In here, we consider two different scenarios as
	in Table 4.4.11:

Table 4 .

 4 4.11: Resistor values for Scenario 2.

	Load	Resistors	Value	Unit
		R 1
	Commercial

Table 4 .

 4 4.12: Simulation results obtained of Scenario 3.

	Load	Power	Power produced [%]	Power consumed [%]	Electricity cost [euros]
		P ug	51.2%	1.9% sold to the UG	4.623
		P es	1.9%	2% for ES charging	-
	Commercial	P loads	-	92% for load usage	-
		P pv	46.9%	-	-
				Total: 4.1%
		P loss	-	R 1 : 0.19%, R 2 : 2.1%,	-
				R 3 : 1.76%, R 4 : 0.05%
		P ug	44%	13.1% sold to the UG	2.857
		P es	6%	6.7% for ES charging	-
	Domestic	P loads	-	78.1% for load usage	-
		P pv	50%	-	-
				Total: 2.1%
		P loss	-	R 1 : 0.2%, R 2 : 0.9%,	-
				R 3 : 0.6%, R 4 : 0.4%
			Time [h]

Table 4 .

 4 4.13: Simulation results obtained of Scenario 3.

	Load	Power	Power produced [%]	Power consumed [%]	Electricity cost [euros]
		P ug	43%	13.1% sold to the UG	2.881
		P es	6.8%	6.8% for ES charging	-
	R 2 =0	P loads	-	77.5% for load usage	-
		P pv	50.2%	-	-
				Total: 2.6%
		P loss	-	R 1 : 1.09%, R 2 : 0%,	-
				R 3 : 0.98%, R 4 : 0.53%
		P ug	43.7%	13.6% sold to the UG	2.765
		P es	3.7%	3.7% for ES charging	-
	R 4 =0	P loads	-	80% for load usage	-
		P pv	52.6%	-	-
				Total: 2.7%
		P loss	-	R 1 : 0.98%, R 2 : 0.61%,
				R

3 : 1.11%, R 4 : 0% -

MVV acronyms stands for Mannheimer Versorgungs und Verkehrsgesellschaft, which is a German energy company.

MV/LV acronyms are for Medium Voltage and Low Voltage respectively.

The duty cycle describes the time for which a signal is active over one period.

The fill factor is defined by the difference between the maximum generated power and the theoretical expected power of a solar cell.

Connecting node is the point where a source or a load is linked to the central network as in Fig.2.2.1 for the connecting nodes 1 . . . 7.

For the sake of convenience, wherever it is straightforward implied by the text, we discard the time dependence.

Keep in mind that the effort e is the voltage and the flow f is the current in the case of electrical circuits, as already mentioned in section 2.3.1.

The circuit's parameters and the values of the resistors are considered as in section 3.3.3:1/C 1sc = α, 1/C 2sc = β, 1/C 3sc = γ, 1/I 1sc = η, 1/I 2sc = µ, 1/R 1sc = ν, 1/R = ζ.

Connecting node is the point where a source or a load is linked to the central network as in Fig. 4.3.4 for the connecting nodes n : 1, n :

2, n : 3, n : 4.

It is the sum of two position vectors X 1 and X 2 which results from the addition of each vector in X 1 to each vector in X 2 .

It is the difference between two position vectors X 1 and X 2 . More precisely, it describes the difference between each vector in X 1 and each vector in X 2 resulting in another vector X 3 : X 3 + X 2 ⊆ X 1 .

The circuit's parameters and the values of the resistors are considered as in section 3.3.2: 1/C 1sc = α, 1/C

2sc = β, 1/C 3sc = γ, 1/C 1b = δ, 1/C 2b = , 1/I 1sc = η, 1/I 2sc = µ, 1/R 1sc = ν, 1/R 1b = σ, 1/R 2b = ψ).

The circuit's parameters and the values of the resistors are considered as in section

3.3.3: 1/C 1sc = α, 1/C 2sc = β, 1/C 3sc = γ, 1/I 1sc = η, 1/I 2sc = µ, 1/R 1sc = ν, 1/R = ζ.

Acknowledgments

It was funded by the French National Research Agency within the framework of the project ANR-15-CE05-004-02 C 3 µ (Components, Control and Communication). During this important period of my life, I have succeeded my objectives due to the help and support of many people and I take this opportunity to express my gratitude to all of them.

List of Tables

More specifically, the tangent flat output is the set of one-forms ω j , j = 1, 2, 3, 4 (the number of one-forms equals to the number of the control inputs of the ES system) and is written as follows: ω 1 = ηp 1sc dp 1sc + µp 2sc dp 2sc + βq 2sc dq 2sc , (3.3.24a)

z 2 = q 3sc + q 1b , (3.3.25b)

)

In the following , the states and the control inputs of the ES system (3.3.18a-3.3.18g) are rewritten in function of the flat outputs (3.3.25a-3.3.25d) and their derivatives:

), (3.3.26a)

)

)

)

)

)

)

Below, the compact flat representation of the states and inputs is deduced in function of the flat outputs and their derivatives, which will be useful later for the explanation of the simulations:

), (3.3.27a)

q 1sc = Φ 3 (z 1 , ż1 , z 2 , ż2 , z2 , z 3 , ż3 , z3 , z 4 , ż4 ), (3.3.27c) q 2sc = Φ 4 (z 4 ), (3.3.27d)

Below, the compact flat representation, as in Definition 3.2.1, of the states and inputs is written:

p 2sc = Φ 45 (z 2 ), (3.4.31b) q 1sc = Φ 46 (z 1 , ż1 , z 2 , ż2 , z 3 , ż3 , z 4 ), (3.4.31c)

)

)

for the second flat output z 2 = p 1sc p 2sc q 2sc q 3sc , the flat representation of the states and the inputs is written below:

)

)

), (3.4.32c)

)

)

Below, the corresponding compact flat representation of the states and inputs is introduced:

)

)

)

) (3.4.28) we have more convergent to the reference value and less noiseaffected results for v s . Contrariwise, i b signal contains more oscillations. In the duty cycles, though, slight disturbances appear for both cases. In general, both sets of flat outputs can be considered for the flat As a first step, the cost function is expressed in terms of the relations below, as they have been introduced in section 2.5.2.2:

where v es (t) and i es (t), v pv (t) and i pv (t), v ug (t) and i ug (t), v loads (t) and i loads (t) denote the voltage and the current on the connecting nodes of the ES, PV, UG and loads system respectively, as in Fig. 4.3.4. The flat representation of the ES system generated by the algorithm in section 3.3.2.1 will be considered. The variable v es (t) can be written in function of the flat outputs as in (3.3.26k). However, this relation is too complicated to be integrated in the optimization problem and written in function of the B-splines. Therefore, an alternative solution will be presented. Taking into account that the converters considered in the model are ideal as already introduced in section 2.4.2.1, the input power of the converter, P sc in (t), is equal to the output power of the converter, P sc out (t), as in (2.4.18), which leads to:

Moreover, the input voltage v sc in of the Split-Pi is always close to 400 V meaning higher than the expecting output voltage, v sc out , coming from the battery. Therefore, the Split-Pi converter will continually operates in down-conversion towards the positive direction and in up-conversion otherwise, which means that d 1sc is equal to 0 and only d 2sc will change values in the (0, 1) interval (for further explanation see also section 4.3.3). Hence, from (2.4.17):

Additionally, through Ohm's law, a connection among the input and output voltage, v sc in (t) and v sc out (t) and current, i sc in (t) and i sc out (t), of the Split-Pi converter and the input voltage, v es (t), and current, i es (t),

Simulation results of the high level with power losses

In this section, the simulation results of the high level are introduced according to (4.3.39a)-(4.3.39h), including the transmission network with power losses. For the simulation, we follow similarly Table 4.4.2 with the added constraints for this case. The profiles generated at the high level for the CU (Fig. 4.4.4a) and the DU load (Fig. 4.4.4b) are illustrated. Moreover, the constraints of the system are considered as previously according to the RPI set found in section 4.4.1.

Concerning the simulations, in the CU case there is a large difference between the P ower Balancing profiles generated from the network without losses (Fig. 4.4.3a) when compared to the network with losses (Fig. 4.4.4a). This happens because in a commercial environment, the consuming hours are during the day. Furthermore, the UG is unable to charge completely the batteries during the night, when the electricity cost is lower, because of the losses that exist in the transmission line (only 0.8% for battery charging). After, during the day, UG and PV sources collaborate (in total 99% power production) to satisfy the consumers' demand (97% consumed from the total power). This is possible considering that we keep same resistor values in the transmission lines, R = 1 Ω. In the case where the batteries are closer to the renewable sources, then the losses among them are lower and the batteries could be more effectively charged. Also, when the distance among the external grid and the microgrid is higher (meaning higher power loss), then the controller gives priority to the parts that are less affected by the power losses. On the other hand, in DU case (Fig. 4.4.3a) where there is no high demand from 12p.m. to 12a.m., the U G and the PV charge the batteries. In general, for the DU profile, similar profiles are obtained in both cases (Fig. 4 performance (note in particular the variation of the electricity cost). Therefore, we cannot ensure which case from Table 4.4.9 is the most effective in order to obtain the optimal profiles for the optimization problem.

On the other hand with differential flatness, apart from the complexity of the flat outputs' calculation to construct the objective function and constraints, we have significant advantages:

we avoid discretization since we solve an optimization problem in continuous time;

since the simulations generated from differential flatness with B-spline parametrization are generated in continuous-time, the results are not affected by discretization approximations or under-sampling. Moreover, the profile is generated in its entirety for the full simulation horizon (and not piece by piece, as would be the case for a discrete optimization problem like MPC); practically, we observed that the total economic cost and required computational resources are lower when exploiting the flat representation;

Appendix A Detailed explanation of the algorithm for the ES, PV and loads systems

Implementation of the algorithm for the ES system

As already indicated in section 3.3.1, the computation hinges on an iteration process that starts from i = 0 and continues step by step, if necessary. Initially, from the linearized system in (3.3.21a)-(3.3.21c), we will introduce the (nm) × n matrices below:

The corresponding matrices of (A.1.1), P 0,[i] and P 1,[i] , at the first step of the algorithm, where i = 0, are equal to 1 :

Then, the dx coordinates are presented through (3.3.21a), (3.3.21b) and (3.3.21c) for i = 0, as shown below:

The time derivative of u [0] will be:

The matrices P ⊥ 1,[0] and w [1] are presented below:

Consequently, the matrix B [0] , considering also (A.1.10a), is as follows:

At this point, the matrix B [0] has full column rank and replaces the previously mentioned matrix B [0] (A.1.6b). Therefore, the iteration process stops. Next, the pseudo-inverse matrix, Z + [0] , is calculated and is equal to:

Therefore, according to (3.3.17), the tangent flat output of the ES system is: For further validation, an approximation of the derivatives is realized by adding a continuous transfer function block after the numerical differentiation in order to filter the signal's response for both sets of flat outputs and avoid the noise created in the data (Fig. A.1.7). We observe that the actual values of the control inputs, v s , d 1sc , d 1sc , i R1b , approximate closely the reference values. Charge q 3sc error 1.5 A.2 PV and loads system A.2.1 Implementation of the algorithm for the Split-pi converter

In this section, we follow the same steps as in Appendix A.1.1 introducing below the corresponding matrices of the algorithm for the Split-Pi converter 2 :

1. The matrices P 0,[0] and P 1,[0] are introduced below:

Next, the dx coordinates are introduced as in (3.3.34) for i = 0, as follows:

with their time derivatives in matrix u [0] presented below:

8. As a consequence, the matrix B [0] , considering also (A.1.11a), is as follows:

which has full column rank, hence the iteration process within the algorithm stops. 9. Next, the pseudo-inverse matrix, Z + [0] , is calculated and is equal to:

(A.2.10) 10. Therefore, according to (3.3.17), the tangent flat output of the Split-Pi converter is:

)

) For further validation, the derivatives are approximated by adding a continuous transfer function block after the numerical differentiation as a filter for both sets of flat outputs in order to mitigate the noise created in the data (Fig.

The previous multiplication in (B.1.1) concludes in a 4-term relation and we proceed to the calculation of each term separately:

Since the T erm 1 results in scalar values, the obtained objective function will be in quadratic form. Similarly, we continue with the calculation of the other terms: where κd + 1 ≤ i, j ≤ κ.