Next Article in Journal
Biogas Produced by Anaerobic Digestion Process and Biodiesel from Date Seeds
Next Article in Special Issue
Stability Analysis of Grid-Forming MMC-HVDC Transmission Connected to Legacy Power Systems
Previous Article in Journal
A Hybrid Simulation Model to Predict the Cooling Energy Consumption for Residential Housing in Hong Kong
Previous Article in Special Issue
Investigating the Converter-Driven Stability of an Offshore HVDC System
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Hierarchical Control Approach for Power Loss Minimization and Optimal Power Flow within a Meshed DC Microgrid

Grenoble INP, Graduate Schools of Engineering and Management, University Grenoble Alpes, LCIS, F-26000 Valence, France
*
Author to whom correspondence should be addressed.
Energies 2021, 14(16), 4846; https://doi.org/10.3390/en14164846
Submission received: 20 June 2021 / Revised: 21 July 2021 / Accepted: 2 August 2021 / Published: 9 August 2021
(This article belongs to the Special Issue Hybrid AC/DC Grid)

Abstract

:
This work considers the DC part of a hybrid AC/DC microgrid with a meshed topology. We address cost minimization, battery scheduling and the power loss minimization within the power distribution network through constrained optimization. The novelty comes from applying differential flatness properties to the microgrid components and formulating the cost and constraints in terms of the associated B-splines parametrization of the flat outputs (the voltages and currents of the system). This allows us to obtain optimal power profiles to minimize the power dissipation and the cost of the electricity purchase from the external grid. These profiles are tracked by a model predictive controller at the higher level, while at a a lower level a controller deals with the operation of the switches within the DC/DC converters. Extensive simulations under nominal and fault-affected scenarios using realistic data validate the proposed approach.

1. Introduction

Microgrids are complex power systems and their efficient operation demands the study on multiple interlocking factors [1], such as power losses, electricity cost, renewable resources, component size and the like. To date, the state of the art in the microgrid topic has shown promising results in the power flow amelioration and performance enhancement [2]. Optimal power flow accounts for the effective and reliable functioning of the system by either minimizing the electricity cost or the energy dissipation. The existence of power loss in the transmission lines can greatly affect the quality of the power distribution. Therefore, the minimization of the power loss constitute a significant factor and leads towards the improvement of the power transmission [3].
Several techniques have been considered recently for the minimization of the power dissipation, emphasizing either the components linked to the microgrid or the transmission-line network [4]. For instance,
Refs. [5,6] focus on the the energy storage (ES) scheduling in combination with the use of renewable sources and their topology within the network. In other works, the researchers seek to reduce the existing power losses which exist in the transmission lines and cannot be avoided or ignored [7]. In the literature, researchers have proposed several methods for optimal power distribution including power loss minimization. In [8], a hierarchical control strategy divided into three levels is proposed to handle the energy management in stand-alone 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 [9] 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. Ref. [10] introduces an optimal power flow algorithm for islanded microgrids to cope with central power dissipation under a centralized three-layer supervision controller.
In the literature and, in particular, in the references already mentioned, the dynamical models are usually taken into account as a group of differential equations without thoroughly representing and studying the power conservation among the components of the system. Regarding the optimal profiles, optimization control methods under constraints are employed using MPC to generate online optimal profiles for discrete time systems penalizing simultaneously the dissipation, the cost and the error that exist among the real and reference profiles. To the best of our knowledge only few works [4,5,6,7,8,9,10,11] deal specifically with the power dissipation issues. The strength of our work is related to the fact that it handles in a multi-objective optimization-based approach both optimal energy management and power losses minimization while satisfying constraints, rejecting disturbances and grasping the system’s multi-timescale characteristics.
Contributions: In [12], the complete three-layer hierarchical control strategy was presented. A constrained optimization-based problem was implemented in the absence of the power losses in the transmission-line network. In this paper, we further develop the optimization problem, particularly of the high level, already presented in [12], taking into account the power dissipation. In the following, we summarize the contributions of this work:
  • Reduction in the power losses existing in the transmission lines is further added to the optimization to enhance the power flow while considering at the same time the cost reduction in the power consumption. The formulation of the high level optimization problem is done through the analytical computation of the system’s flat representation together with B-spline parametrization;
  • Thorough presentation of the central transmission network describing in detail its dynamical model. The power dissipation is considered through the voltage drops among the four connecting nodes n: 1 , 2 , 3 , 4 [13], illustrated in Figure 1 and the corresponding transmission lines are represented through four resistances, R 1 , 2 , 3 , 4 . Additional constraints are taken into consideration to preserve the DC voltage close to the reference equal to 400 V;
  • Verification of the meshed topology of the central network in the case of unexpected events (e.g., blackouts).
Note that the objective function will be posed in a flat representation (rewritten in terms of a flat output and its derivatives). Likewise will be done for the constraints. Further parametrizing the flat output by a family of B-spline functions will allow to check continuous time constraints. More details on the concepts of flatness and B-spline parametrization can be found in [14,15,16].
Outline:Section 2 concerns the abbreviations. In Section 3 a brief description of the ES (Energy Storage) dynamics is introduced and the detailed model of the central transmission network are introduced. Section 4 contains the multilevel hierarchical control problem, with the high and the middle levels being detailed in particular. The results are introduced in Section 5. Finally, Section 6 presents the outcomes of this work.

2. Notations

E S , e s energy storage system
U G , u g utility grid
P V , p v solar panel system
D E R s distributed energy resources
x ( t ) state vector of the system R n
u ( t ) input vector of the system R m
y ( t ) output vector of the system R m
nconnecting node between a source or a load and the DC-bus
K i B a M Kinetic Battery Model
Pelectrical power
S w switch of the DC/DC converter
d 1 s c , d 2 s c duty cycles of the switches within the Split Pi converter
q 1 b , q 2 b available and bound charge of the KiBaM battery
bKiBaM battery
s c Split−Pi Converter
Ccapacitor
Iinductor
R 1 b resistance between the KiBaM battery and the Split−Pi
. converter
R 2 b resistance of the KiBaM battery
R 1 , . R 2 , . R 3 , . R 4 resistances of the transmission network
R 1 s c resistance among the DC network and the Split−Pi
. converters
i s c _ o u t output current of the Split−Pi converter
v s c _ o u t output voltage of the Split−Pi converter
i s c _ i n input current of the Split−Pi converter
v s c _ i n input voltage of the Split−Pi converter
rderivative of the B-spline
p i i th control point
b i , d 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
S κ , d 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 functions
Tknot vector R N + d
τ κ κ th knot
κ number of knots
eelectricity price
Q y ˜ , . R u ˜ weight matrices

3. DC Microgrid Model Description

The DC microgrid, part of a hybrid AC/DC architecture with a meshed topology (Figure 1), considered in this work is also described in [12] and is composed by a set of PVs, an ES composed by lead acid batteries and a set of loads. These components are linked to the central transmission network (around 400 V DC) through a particular type of DC/DC converters: Split-Pi converters, which regulate the input/output voltage from and towards the sources and the loads. The DC microgrid is connected to the UG through an AC/DC converter. In this work, we analyze the dynamical model of the DC-bus together with the ES system (Split-Pi converter connected to the battery). However, the following approach can be similarly applied for the primarily proposed design in [17]. The overall hybrid AC/DC grid contains an AC/DC converter, placed between the grid and the DC microgrid (Figure 1), where supplementary buses (AC or DC) or sources (AC or DC) could be included, since a meshed topology of a multiple-line transmission network is considered. However, for an AC/DC converter, an AC bus or an AC source, the characteristics of the AC power should be taken into account, such as the active power, the reactive power, the frequency and the like [18].
This section builds upon our previous paper [12], where the mathematical model of the ES system was explicitly introduced. In this work, we will continue with the explicit description of the dynamics of the central transmission network, since the power losses are now considered and resistors exist to replace every transmission line. The mathematical methodology employed relies upon the PH formalism [19]. This method provides an explicit description of the routing elements (inductors, capacitors, resistors) existing in the electrical network. An advantage of this method is that the set of equations can be directly derived by the associated Bond Graph of the electrical system [20,21], which is a graph-oriented modeling tool for the dynamics of convoluted physical systems. The modeling methodology are thoroughly presented and analyzed in [12,17].

3.1. Brief Description of the ES Subsystem

The ES constitutes the Split-Pi converter connected to the battery as in Figure 2. The circuit of the battery is structured according to the Kinetic Battery Model (KiBaM), consisting of two capacitors and a resistance between them. The charges of the capacitors represent the total charge of the battery. Concerning the Split-Pi converter, it is a bidirectional buck-boost converter able to regulate the voltage in both directions.
The state-space representation of the ES system was already explicitly presented in [12] concluding to the following nonlinear state-space mathematical model:
p ˙ 1 s c ( t ) = q 1 s c ( t ) C 1 s c q 2 s c ( t ) C 2 s c ( 1 d 1 s c ( t ) ) ,
p ˙ 2 s c ( t ) = q 2 s c ( t ) C 2 s c ( 1 d 2 s c ( t ) ) q 3 s c ( t ) C 3 s c ,
q ˙ 1 s c ( t ) = v D C R 1 s c q 1 s c ( t ) C 2 s c R 1 s c p 1 s c ( t ) I 1 s c ,
q ˙ 2 s c ( t ) = p 1 s c ( t ) I 1 s c ( 1 d 1 s c ( t ) ) p 2 s c ( t ) I 2 s c ( 1 d 2 s c ( t ) ) ,
q ˙ 3 s c ( t ) = p 2 s c ( t ) I 2 s c i R 1 b ( t ) ,
q ˙ 1 b ( t ) = i R 1 b ( t ) q 1 b ( t ) C 1 b R 2 b + q 2 b ( t ) C 2 b R 2 b ,
q ˙ 2 b ( t ) = q 1 b ( t ) C 1 b R 2 b q 2 b ( t ) C 2 b R 2 b .
The p and q variables correspond to the magnetic flux of the inductors and the charges of the capacitors, respectively. These variables ( p 1 s c , 2 s c , q 1 s c , 2 s c , 3 s c ) structure the state vector of the ES system, x e s R 7 × 1 . Next, v D C and i R 1 b form the input vector u e s R 2 × 1 of the ES and i D C and v R 1 b form the output vector y e s R 1 × 2 . The circuit’s parameters are I 1 s c , I 2 s c , C 1 s c , C 2 s c , C 3 s c , C 1 b , C 2 b , with I denoting the inductance of the inductor and C the capacitance of the capacitor. The state vector, the input vector and the output vector of the system are given analytically in [12]. Furthermore, from the ES circuit in Figure 2 and [12], we have the following relation which will be useful later in the energy management:
i b ( t ) = i R 1 b ( t ) ,
v b ( t ) = q 1 b ( t ) C 1 b .
The Split-Pi converter is considered as an ideal one ( see also [13]), that is, the following relation is satisfied: P s c _ i n = P s c _ o u t . Ohm’s law is taken into account to link the input and output voltage and current of the Split-Pi converter as in [12,13]. Hence, the ES power is given as in [13]: P e s ( t ) = v e s ( t ) i e s ( t ) = ( v b ( t ) + i b ( t ) R 1 b ) i b ( t ) + R 1 s c ( α ( t ) i b ( t ) ) 2 and will be used in the analysis and calculation of the cost function.

3.2. Dynamical Representation of the DC Bus

The dynamical model of the four-line transmission network (Figure 1) is developed by its Bond graph as already explained in [17]. Ref. [17] provides the Bond graph of the central transmission network taking into account each transmission lines as an RL circuit. In this work, the central network is further simplified eliminating the inductors. In DC networks, the existence of the inductors is necessary only when the current flow of the system varies. However, once the system comes to its steady state, where the current flow remains stable, the inductors have no effect in the DC networks. Hence, each transmission line will be considered as a resistor and the power conservation equation is as follows (Figure 1):
P p v ( t ) P e s ( t ) + P u g ( t ) P l o a d s ( t ) P R ( t ) = 0 .
P p v , e s , u g , l o a d s is the power of the sources and the loads, while P R is equal to P R 1 + P R 2 + P R 3 + P R 4 the power loss for the transmission lines R 1 , R 2 , R 3 , R 4 . The correspondingly-modified Bond graph of the central DC bus is depicted in Figure 3.
Hence, from the the Bond graph in Figure 3, where only dissipative elements are included, the following is obtained:
i p v i u g i e s i l o a d s = 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 p v v u g v e s v l o a d s .
Therefore, from the above mentioned equations, the current flows on the connecting nodes are deduced as in [13]: (i) n : 1 i l o a d s = i R 1 + i R 2 = v e s v l o a d s R 1 + v u g v l o a d s R 2 ; (ii) n : 2 i e s = i R 3 i R 1 = v p v v e s R 3 v e s v l o a d s R 1 ; (iii) n : 3 i u g = i R 4 i R 2 = v p v v u g R 4 v u g v l o a d s R 2 ; (iv) n : 4 i p v = i R 3 i R 4 = v p v ( t ) v e s ( t ) R 3 v u g ( t ) v p v ( t ) R 4 . The goal is to express all the powers appearing in the grid in terms of voltages and resistances (i.e., to avoid the explicit appearance of currents). Through the relations of the current flows on the connecting nodes previously presented we make use of the fact that the power in every electrical element is P = i × v . Hence, every power variable of the DC network is expressed in function of the voltages on the connecting nodes (see also Figure 1) as in [13]: (i) P p v = v p v i p v = v p v [ i R 3 i R 4 ] = v p v v p v v e s R 3 v u g v p v R 4 ; (ii) P e s = v e s i e s = v e s [ i R 3 i R 1 ] = v e s v p v v e s R 3 v e s v l o a d s R 1 ; (iii) P u g = v u g i u g = v u g [ i R 4 i R 2 ] = v u g v u g v p v R 4 v l o a d s v u g R 2 ; (iv) P l o a d s = v l o a d s i l o a d s = v l o a d s [ i R 1 i R 2 ] = v l o a d s v e s v l o a d s R 1 v l o a d s v u g R 2 . For the power dissipation within the central transmission network, P R 1 , P R 2 , P R 3 , P R 4 , we use the other equivalent power formulation, P = v 2 × R , which allows deduction of the following relations:
P R 1 = [ v e s v l o a d s ] 2 R 1 ,
P R 3 = [ v p v v e s ] 2 R 3 ,
P R 2 = [ v l o a d s v u g ] 2 R 2 ,
P R 4 = [ v u g v p v ] 2 R 4 .
The aforementioned equations are significant and will be proved useful later in structuring the objective function in order to handle the power loss minimization problem.

4. Optimization Objectives and Constraints

Hereinafter, the objectives and constraints of the energy management optimization problem will be introduced.

4.1. Objectives

The main objective of this work is the optimal power flow while minimizing the power dissipation in the DC bus (Figure 4). This problem aims to optimize the power flows within the microgrid so as to minimize the power losses during transmission with the minimum electricity cost. Therefore, the general cost function which penalizes the energy dissipation will be divided in two parts:
  • 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 0 t f e ( t ) P u g = t 0 t f e ( t ) ( P l o a d s + P e s P p v ) ,
    where reference profiles will be taken into account for the PV and the loads. The P u g is replaced, exploiting the power conservation law ( P u g = P l o a d s + P e s P p v ), without considering the power loss (3);
  • The optimal power flow, minimizing the power dissipation in the DC bus. Therefore, the cost function will be the following:
    J 2 ( d ( t ) ) = t 0 t f P R .
    where P R corresponds to the power dissipation inside the central network (5)–(8). Relation (3) 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 here and in [12].

4.2. Constraints

At this point, the general constraints considered at the energy management problem will be analyzed. The ES system constitutes a major component of the overall scheme. Concerning the batteries, their limited lifetime demand to respect certain constraints. Hence, the constraints for the ES system are (see also Figure 2):
q 1 b m i n q 1 b ( t ) q 1 b m a x , q 2 b m i n q 2 b ( t ) q 2 b m a x , Charging   mode : v b , c h a r g i n g m i n v b , c h a r g i n g ( t ) v b , c h a r g i n g m a x , i b , c h a r g i n g m i n i b , c h a r g i n g ( t ) i b , c h a r g i n g m a x , Dischaging   mode : v b , d i c h a r g i n g m i n v b , d i c h a r g i n g ( t ) v b , d i c h a r g i n g m a x , i b , d i c h a r g i n g m i n i b , d i c h a r g i n g ( t ) i b , d i c h a r g i n g m a x .
Furthermore, the voltage on the connecting nodes should always be close to v D C = 400 V. Therefore, the constraints of the connecting nodes are mentioned as follows ( v D C m i n , h and v D C m a x , h will be the limits of the DC voltage on every connecting node which will be close to the nominal, v D C = 400 V) (Figure 1):
v D C m i n , h v u g ( t ) v D C m a x , h ,
v D C m i n , h v p v ( t ) v D C m a x , h ,
v D C m i n , h v e s ( t ) v D C m a x , h ,
v D C m i n , h v l o a d s ( t ) v D C m a x , h .
Additionally, constraints for the duty cycles d 1 s c and d 2 s c of the Split-Pi converter are considered as shown below:
0 d 1 s c 1 ,
0 d 2 s c 1 .
Finally, maximum and minimum limits for the external grid power generation must be considered:
P u g m i n P u g ( t ) P u g m a x ,
with P u g ( t ) = v u g ( t ) i u g ( t ) .

ES System Flat Representation

In this subsection, for the ES flat representation (1a)–(1g), the flat outputs are already described in our previous work in [12] and are recalled here for a clear understanding of the proposed approach: (i) z 1 ( t ) = 1 I 1 s c p 1 s c ( t ) 2 2 + 1 I 2 s c p 2 s c ( t ) 2 2 + 1 C 2 s c q 2 s c ( t ) 2 2 ; (ii) z 2 ( t ) = q 3 s c ( t ) + q 1 b ( t ) ; (iii) z 3 ( t ) = q 2 b ( t ) ; (iv) z 4 ( t ) = q 2 s c ( t ) . Replacing the aforementioned flat outputs into the PH model (1a)–(1g), the flat representation of the system is acquired, as well as the flat representation of the states and the inputs of the system.
Since the ES system has been already written in terms of the flat outputs, the ES voltage and the current, v e s ( t ) and i e s ( t ) , can be written in function of a family of N B-splines functions of order d which are denoted as b i , d ( t ) . These basis functions are summed together (and weighted through their control points p i ) in order to describe the flat output z 3 ( t ) . Hence, writing the flat output z 3 in function of the B-spline curves, z 3 ( t ) = i = 1 N p i b i , d ( t ) = P B d ( t ) , we can now express (15a) and (15b) by exploiting the derivation property of B-splines (which allows characterization of derived splines in terms of lower order splines; matrices M d , d 1 , S k , d 1 , d are computed accordingly as may be seen for example in [22]). Through the B-spline derivation property, we parametrize firstly the voltage of the battery, v b , and the current of the battery, i b , in function of which the ES voltage and the current, v e s ( t ) and i e s ( t ) , are written: (i) v e s ( t ) = v b ( t ) + i b ( t ) R 1 b α ( t ) + R 1 s c α ( t ) i b ( t ) ; (ii) i e s ( t ) = α ( t ) i b ( t ) [13]. Therefore, we have the following:
p i v b = 1 C 2 b p i + R 2 b P M d , d 1 S κ , d 1 , d i ,
p i i b = 1 + C 1 b C 2 b P M d , d 1 S κ , d 1 , d i + C 1 b R 2 b · P M d , d 2 S κ , d 2 , d i ,
where p i i b and p i v b are the control points for the battery’s current and voltage, respectively. Next, for v b ( t ) and i b ( t ) as already presented in [12], we have:
v b ( t ) = b i , d ( t ) i = 1 N 1 C 2 b p i + R 2 b P M d , d 1 S κ , d 1 , d i ,
i b ( t ) = b i , d ( t ) i = 1 N 1 + C 1 b C 2 b P M d , d 1 S κ , d 1 , d i + C 1 b R 2 b P M d , d 2 S κ , d 2 , d i , a t τ κ τ κ + 1 .
All the variables and the parameters of the Equations (16a) and (16b) are introduced and explained in the notations section.

5. Hierarchical Control Approach with Constrained Optimization

This section presents the hierarchical control problem analytically, describing the high and the middle level including the power dissipation. The high level profiles for the ES, P e s , and the UG, P u g , are generated by mitigating the electricity cost and the power losses. At the high and the middle level, the whole dynamics of the ES system is taken into account, as will be shown later, in order to link the power losses of the central transmission network to the battery dynamics. Afterwards, the acquired high level profiles are used as references at the middle level for tracking under perturbation. The discretized model of the battery is considered, already presented in [12], and the main objective is to minimize the deviations amongst the reference and the real 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 e s , are applied to the ES PH model to control the switches.
As aforementioned, the flatness-based optimization approach of [12] is reformulated adding the power losses within the central network as in Figure 1. Therefore, the dynamics of the system change 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. In the present case, since the power losses in the central transmission lines are considered, the voltage in the common DC-bus will not be stable and it will vary regarding the constraints of (11a)–(11d). On each connecting node 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 Figure 5, the control scheme of the energy management problem, under the presence of the power dissipation, is illustrated. Optimal profiles are generated for the ES system, the duty cycles through the α factor and the voltages on the connecting nodes.

5.1. High Level

This section will introduce the optimization problem of the high level including details for its flat representation and B-spline parametrization for continuous-time constraint validation.

Definition of the Optimization Problem with the Power Loss Included

The principal objective of the high level problem is the generation of reference trajectories in continuous time, exploiting the flat output representation of the ES system. 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 0 t f e ( t ) P u g d t ,
where e ( t ) is the electricity price.
At first, the power losses are taken into consideration (3), hence P u g ( t ) is equal to:
P u g ( t ) = P l o a d s ( t ) + P R ( t ) + P e s ( t ) P p v ( t ) P e s ( t ) = P l o a d s ( t ) P R ( t ) + P p v ( t ) + P u g ( t ) ,
Hereinafter, the objective function is analytically calculated and has the following form:
J = t 0 t f e ( t ) P l o a d s ( t ) + P R ( t ) + P e s ( t ) P p v ( t ) d t
and the complete optimization problem becomes:
m i n α , v u g , v p v , v e s , v b , i b t 0 t f e ( t ) [ P l o a d s ( t ) + P R ( t ) P p v ( t ) + P e s ( t ) i b ( t ) v b ( t ) ( t ) ] d t ,
s u b j e c t a t o : . t h e . s y s t e m . d y n a m i c s . ( 1 a ) ( 1 g ) ,
. t h e . p o w e r . c o n s e r v a t i o n . ( 18 ) ,
v D C m i n , h v l o a d s ( t ) , v e s ( t ) , v p v ( t ) v u g ( t ) v D C m a x , h ,
v b m i n , h v b ( t ) v b m a x , h ,
i b m i n , h i b ( t ) i b m a x , h ,
q 2 b m i n , h q 2 b ( t ) q 2 b m a x , h ,
P u g m i n , h + P p v ( t ) P l o a d s ( t ) P R ( t ) P e s ( t ) P e s ( t ) P u g m a x , h P p v ( t ) + P l o a d s ( t ) P R ( t ) .
Additionally, the following are also considered to write the final form of the objective function of (20a):
  • the voltage, v l o a d s , on the connecting node 1 is written in function of the input voltage, v e s , and the input current, i e s , of the ES system:
    v l o a d s ( t ) = R 1 i e s ( t ) + ( 1 + R 1 R 3 ) v e s ( t ) R 1 R 3 v p v ( t ) ;
  • the constraint for the consumer’s demand is considered as follows:
    P l o a d s ( t ) ϵ l o a d s v l o a d s ( t ) + v e s ( t ) R 1 v u g ( t ) + v l o a d s ( t ) R 2 v l o a d s ( t ) P l o a d s ( t ) + ϵ l o a d s ,
    where ϵ l o a d s is the soft constraint to relax the consumers’ demand in order to guarantee that the optimization problem will be resolvable.
The voltages on the connecting nodes, v e s , v u g , v l o a d s and v p v , obey to (11a)–(11d). Finally, the α ( t ) factor and the P u g ( t ) (taking into account also (18)) will be:
0 < α ( t ) < 1 ,
P u g m i n , h P u g ( t ) P u g m a x , h
P u g m i n , h P l o a d s ( t ) + P e s ( t ) P p v ( t ) + P R ( t ) P u g m a x , h .

5.2. Flat Representation and B-Spline Parametrization of the Optimization Problem

Next, the optimization problem in (20a)–(20h) is revised in terms of the flat outputs. The cost function is expressed in terms of the relations in (5)–(8). The flat representations of the v e s ( t ) and i e s ( t ) variables of the ES system are deduced in terms of the flat representations of the battery’s current and voltage, i b ( t ) and v b ( t ) , as previously presented in [12]. Substituting (5)–(8) in (19), the cost function becomes:
J = t 0 t f e ( t ) [ Q c o s t P l o a d s ( t ) + P e s ( t ) v e s ( t ) i e s ( t ) P p v ( t ) + Q l o s s ( v l o a d s ( t ) + ( v e s ( t ) ) 2 R 1 + + v u g ( t ) + ( v l o a d s ( t ) ) 2 R 2 + v p v ( t ) + ( v u g ( t ) ) 2 R 4 + v e s ( t ) + ( v p v ( t ) ) 2 R 3 ) ] d t .
Since i e s ( t ) and v e s ( t ) are in function of the battery’s current and voltage, i b and v b , they can be parametrized in terms of the B-splines according to (16a) and (16b). Similarly, (23a) and (23b) must be also parametrized in terms of the B-splines. The α factor is, firstly, parametrized considering 1 α ( t ) ( 1 , + ) and a group of B-splines basis functions with d α order as it is written below:
α ( t ) 1 = j = 1 N α p j α b i , d α ( t ) .
The number of control points p i α is deduced as N α . Then, applying Theorem A1 from Appendix A, a sufficient condition for ensuring (20a) is that each of the scalar control points are supraunitary:
p i α > 1 , i = 1 , , N α .
Since α lies in the interval of ( 0 , 1 ) , meaning that it is positive, the following can be considered for the second part of v e s :
R 1 s c | i b ( t ) | R 1 s c α ( t ) i b ( t ) R 1 s c | i b ( t ) | .
Therefore, the constraint v D C m i n , h v e s ( t ) v D C m a x , h , where v e s is considered, is valid if and only if:
v b ( t ) + i b ( t ) R 1 b α ( t ) R 1 s c | i b ( t ) | v D C m i n , h ,
v b ( t ) + i b ( t ) R 1 b α ( t ) + R 1 s c | i b ( t ) | v D C m a x , h .
Then, for (11a)–(11d) and (23b)), v e s , i e s and P e s are combined to arrive at the following:
p i v b + R 1 b p i i b p j α R 1 s c | p i i b | v D C m i n , h ,
p i v b + R 1 b p i i b p j α + R 1 s c | p i i b | v D C m a x , h ,
where p κ , i v b and p κ , i i b are defined by (15a) and (15b) and κ , i and j N satisfy d 1 κ n 1 , κ d + 2 i κ + 1 and 1 j N α . After various computations, which is given explicitly in Appendix B, the previous inequalities are replaced by the sufficient formulations:
v D C m i n , h i = 1 N j = 1 N α ( p i v b + R 1 b p i i b ) p j α b i , j , d ( t ) i = 1 N j = 1 N α R 1 s c | p i i b | p j α b i , j , d ( t ) i = 1 N j = 1 N α p i v b + R 1 b p i i b R 1 s c | p i i b | p j α b i , j , d ( t ) ,
v D C m a x , h i = 1 N j = 1 N α ( p i v b + R 1 b p i i b ) p j α b i , j , d ( t ) + i = 1 N j = 1 N α R 1 s c | p i i b | p j α b i , j , d ( t ) i = 1 N j = 1 N α ( p i v b + R 1 b p i i b ) + R 1 s c | p i i b | p j α b i , j , d ( t ) .
Concerning the P u g constraint and taking into account the i u g current flow on the connecting node n : 1 as previously mentioned, it is defined below:
P u g ( t ) = i u g ( t ) v u g ( t ) = v u g ( t ) [ i R 4 ( t ) i R 2 ( t ) ] = v u g ( t ) ( v p v ( t ) + v u g ( t ) R 4 v u g ( t ) + v l o a d s ( t ) R 2 ) = v u g 2 ( t ) 1 R 2 + 1 R 4 v l o a d s ( t ) R 2 + v p v ( t ) R 4 v u g ( t ) .
Replacing v l o a d s ( t ) with (21) leads to:
P u g ( t ) = v u g 2 ( t ) 1 R 4 + 1 R 2 1 R 4 v u g ( t ) v p v ( t ) v u g ( t ) [ R 1 i e s ( t ) + + 1 + R 1 R 3 v e s ( t ) R 1 R 3 v p v ( t ) ] .
Including also the relations amongst v e s and i e s with v b and i b , as presented also in [13], P u g ( t ) is denoted as:
P u g ( t ) = 1 R 2 + 1 R 4 v u g 2 ( t ) + R 1 R 3 R 2 1 R 4 v p v ( t ) v u g ( t ) v u g ( t ) 1 R 2 α ( t ) i b ( t ) R 1 + 1 + R 1 R 3 v b ( t ) + R 1 b i b ( t ) α ( t ) + α ( t ) i b ( t ) R 1 s c = = v u g 2 ( t ) 1 R 2 + 1 R 4 + v p v ( t ) v u g ( t ) R 1 R 2 R 3 1 R 4 v u g ( t ) 1 R 2 R 1 R 3 + 1 v b ( t ) + i b ( t ) R 1 b α ( t ) + R 1 R 1 s c R 3 + R 1 + R 1 s c α ( t ) i b ( t ) .
Therefore, finally, the constraint (23b) is defined as:
P u g m a x , h ( 1 R 2 + 1 R 4 ) v u g 2 ( t ) + v p v ( t ) v u g ( t ) ( 1 R 4 + R 1 R 3 R 2 ) v u g ( t ) 1 R 2 ( 1 + R 1 R 3 ) v b ( t ) + i b ( t ) R 1 b α ( t ) | i b ( t ) | ( R 1 + R 1 s c + R 1 s c R 1 R 3 ) ,
P u g m i n , h ( 1 R 2 + 1 R 4 ) v u g 2 ( t ) + v p v ( t ) v u g ( t ) ( R 1 R 3 R 2 1 R 4 ) 1 R 2 v u g ( t ) ( 1 + R 1 R 3 ) v b ( t ) + i b ( t ) R 1 b α ( t ) + | i b ( t ) | ( R 1 s c R 1 + R 1 s c + R 1 R 3 ) .
Additionally, considering that | i b ( t ) | α ( t ) i b ( t ) | i b ( t ) | and v b ( t ) + R 1 b i b ( t ) α ( t ) > 0 we can rewrite the two previous constraints, (34a) and (34b), as follows:
( 1 + R 1 R 3 ) v b ( t ) + i b ( t ) R 1 b α ( t ) | i b ( t ) | ( R 1 s c + R 1 + R 1 s c R 1 R 3 ) R 2 v u g ( t ) [ P u g m a x , h v u g 2 ( t ) ( 1 R 2 + 1 R 4 ) v p v ( t ) v u g ( t ) ( R 1 R 3 R 2 1 R 4 ) ] ,
( 1 + R 1 R 3 ) v b ( t ) + i b ( t ) R 1 b α ( t ) + | i b ( t ) | ( R 1 s c + R 1 + R 1 s c R 1 R 3 ) R 2 v u g ( t ) [ P u g m i n , h v u g 2 ( t ) ( 1 R 2 + 1 R 4 ) v p v ( t ) v u g ( t ) ( R 1 R 3 R 2 1 R 4 ) ] .
Next, through (A1a) and (A1b) in Appendix B, the left part of (35a) and (35b) is defined in terms of the B-splines:
( R 1 R 3 + 1 ) v b ( t ) + i b ( t ) R 1 b α ( t ) ± | i b ( t ) | ( R 1 s c + R 1 + R 1 s c R 1 R 3 ) = ( 1 + R 1 R 3 ) i = 1 N j = 1 N α ( p i v b + R 1 b p i i b ) p j α b i , j , d ( t ) ± ( R 1 + R 1 s c + R 1 R 1 s c R 3 ) i = 1 N j = 1 N α | p i i b | b i , j , d ( t ) ,
which is proven similarly as in (A2a) and (A2b).
Hence, the optimization problem of (20a)–(20h) becomes:
m i n p , p α , v p v ( t ) , v u g ( t ) t 0 t f e ( t ) [ Q c o s t [ ( v b ( t ) + R 1 b ) i b ( t ) 2 + R 1 s c ( i b ( t ) ) 2 α ( t ) P p v ( t ) +
+ P l o a d s ( t ) ] + Q l o s s [ ( i b ( t ) R 1 b + v b ( t ) α ( t ) + v l o a d s ( t ) ) 2 + R 1 s c α ( t ) i b ( t ) R 1 +
+ ( v u g ( t ) + v l o a d s ( t ) ) 2 R 2 + ( v p v ( t ) i b ( t ) R 1 b + v b ( t ) α ( t ) + R 1 s c i b ( t ) ) 2 α ( t ) R 3 +
+ ( v u g ( t ) v p v ( t ) ) 2 R 4 ] ] d t
s u b j e c t a t o : . t h e . s y s t e m . d y n a m i c s . ( 1 a ) ( 1 g ) ,
. t h e . p o w e r . c o n s e r v a t i o n . ( 18 ) ,
t h e a v o l t a g e a c o n s t r a i n t s a o n a t h e a c o n n e c t i n g a n o d e s a ( 11 )
v b m i n , h i = 1 N b i , d ( t ) p i v b v b m a x , h ,
i b m i n , h i = 1 N b i , d ( t ) p i i b i b m a x , h ,
q 2 b m i n , h i = 1 N b i , d ( t ) p i q 2 b m a x , h ,
t h e a P u g a c o n s t r a i n t s a ( 34 a ) , ( 34 b ) ,
t h e a c e n t r a l a n e t w o r k a r e l a t i o n s a c o n s t r a i n t s a ( 21 ) , ( 22 ) .
i b ( t ) and v b ( t ) are deduced in function of the B-splines as in (16a) and (16b), respectively. Q c o s t and Q l o s s are the weight matrices corresponding to the cost minimization and power loss mitigation.

5.3. Middle Level

Next, the middle level of the hierarchical problem will be described. The battery current, i b , battery voltage, v b , and input voltage of the ES system v e s reference profiles from the high level will be considered here as i b r e f , v b r e f and v e s r e f correspondingly. Likewise for the α factor which will be referred to as α r e f at the middle level.
Therefore, we introduce a tube-MPC controller to follow the Split-Pi output voltage reference profile, v s c _ o u t r e f , under perturbation. The output voltage reference of the converter is written in function of the battery current and voltage reference profiles from the high level [13]:
v s c _ o u t r e f ( t ) = v b r e f ( t ) + i b r e f ( t ) R 1 b ,
according to Ohm’s law (Figure 2). We proceed to the battery’s discretization via the standard Euler discretized as already introduced in [12].
At this point, we continue using the tube-MPC controller and the discretized dynamics of the system. However, since the power losses are taken into account, the MPC tracking problem is reformulated as follows:
min u ˜ ( k ) i = k k + N p 1 ( v ˜ e s ( i ) v ˜ e s r e f ( i ) ) Q v ˜ e s ( v ˜ e s ( i ) ) + + ( u ˜ ( i ) u ˜ r e f ( i ) ) R u ˜ u ˜ ( i ) u ˜ r e f ( i )
s u b j e c t a t o : . t h e . s y s t e m . d i s c r e t i z e d . d y n a m i c s ,
v ˜ b m i n , m v ˜ b ( k ) v ˜ b m a x , m ,
i ˜ b m i n , m i ˜ b ( k ) i ˜ b m a x , m ,
q ˜ 2 b m i n , m q ˜ 2 b ( k ) q ˜ 2 b m a x , m ,
v ˜ e s m i n , m v ˜ e s ( k ) v ˜ e s m a x , m ,
P ˜ u g m i n , m P ˜ u g ( k ) P ˜ u g m a x , m .
In the last constraint, P u g ( t ) , the power losses must be considered in the power conservation equation as in (18) in Section 5.1. Therefore, (39g) is replaced by:
P ˜ u g m i n , m + P ˜ p v ( k ) P ˜ l o a d s ( k ) P ˜ R ( k ) P ˜ e s ( k ) ,
P ˜ u g m a x , m + P ˜ p v ( k ) P ˜ l o a d s ( k ) P ˜ R ( k ) P ˜ e s ( k ) ,
P ˜ e s ( k ) is equal to i ˜ e s ( k ) v ˜ e s ( k ) and P ˜ R ( k ) is equal to P ˜ R 1 ( k ) + P ˜ R 2 ( k ) + P ˜ R 3 ( k ) + P ˜ R 4 ( k ) . v ˜ e s ( k ) and i ˜ e s ( k ) are computed regarding the α r e f obtained at the high level as in [13]:
v e s ( k ) = v b ( k ) + i b ( k ) R 1 b α r e f , h ( k ) + R 1 s c α r e f ( k ) ,
i e s ( k ) = α r e f ( k ) i b ( k ) .
At this point, we mention that the low level remains the same as in [12] since it is a local controller concerning the operation of the Split-Pi converter locally and its dynamics are not influenced by the addition of the power dissipation in the DC. The single difference is the input voltage, v e s , of the ES system on the connecting node 2, which varies among the considered constraints (11a)–(11b). The tracking profiles obtained for the voltage, v e s , at the middle level are considered.

6. Simulation Results

In this section, the simulation results of the of the high, the middle and the low levels are presented. The DC microgrid parameters, considered for the simulations, are illustrated in Table 1. The variables and constraints for the high level and the middle level are depicted in Table 2. Reference profiles for the PV and the loads are considered, obtained considering real irradiation and temperature data. Additionally, for the ES system, we take into account a set of lead acid batteries (AGM 12-165). A DC breaker is used to link the DC microgrid with the UG. We use MATLAB 2015 a for the simulations and YALMIP optimization toolbox [23] for both high and middle levels in order to be able to use the IPOPT solver [24] which handles nonlinear optimization problems. Concerning the low level, the PH ES system was designed in MATLAB/Simulink with the goal to confirm the appropriate switching activity within the Split-Pi [12].

6.1. Simulation Results of the High Level

In this section, the simulation results of the high level are introduced according to (37a)–(37l). The profiles generated at the high level for the commercial use (CU) (Figure 6a), where the load is high during the day, and the domestic use (DU) (Figure 6b), where the load is high during the afternoon, are illustrated. Moreover, the constraints of the system are considered according to the RPI set found in [12].
Concerning the simulations, in the CU case for the P o w e r B a l a n c i n g profiles (Figure 6a) the consuming hours are during the day. The UG cannot completely charge the ES at night, when the cost of the electricity cost is low, because of the power loss existing in the transmission lines (only 0.8 % for battery charging). Afterwards, 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, the losses among them are lower and the batteries could be more effectively charged. Moreover, when the distance among the external grid and the microgrid is higher (meaning higher power loss), then the controller gives priority to the parts which are less affected by the power losses. On the other hand, in DU case the PV and the U G charge the ES from 12 p.m. to 12 a.m. when the demand is low. In the afternoon an increase in the demand is observed, while the rest of the day the demand remains broadly stable. As aforementioned, there is, also, enough energy to sell to the UG, close to 13 % of the total consumed power. Hence, to make profit in a commercial environment, the use of larger installations for the renewable resources is important.
Likewise, the batteries’ reference profiles for the voltage, current and state of charge are introduced where the constraints are verified. The electricity cost is also calculated which is 4.318 euros for the CU and 2.713 euros for the DU. The cost without using the ES system is 4.737 euros for the CU profile and 4.732 euros for the DU profile, where we observe that in the first case (CU profile) 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.
Then, Figure 7 depicts the power losses and the constraint validation for the voltage on the connecting nodes. In the CU profile (Figure 7a), 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 (Figure 7b). 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 3 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.
The previous 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.

6.2. Simulation Results of the Middle and Low Level

At this section, the middle level results are introduced considering the optimal profiles generated at the high level for the battery current, i b , and the voltage, v b as references (Figure 6). Table 2 depicts the parameters for the simulations of the middle level. Here, MPC is used for reference tracking considering a prediction horizon N p equal to 5 and a sampling time T s equal to 5 min. The system, at this point, is discretized and it follows the optimization problem under constraints as in (39a)–(39g) with losses in the transmission network. Additionally, a perturbation is added to the system that is always lower than the difference between the maximum and minimum state value defined by the RPI set in [12].
Figure 8 shows the middle level profiles of the energy management problem, and the control input, v s c _ o u t of the Split-Pi converter, in terms of i b and v b as in (38). From the figures and the Table 4, 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.
Furthermore, there are 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 or the DU demand.
Next, the low level results are presented (Figure 9) tracking the profiles obtained at the middle level for the i b and the v b under perturbation (Figure 8). 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 [12]. Apart from the voltage and current reference profiles of the battery, the reference profile of the voltage, v e s , entering from the central transmission network is also considered (see Figure 7a,b).

6.3. Optimal Profile Generation of Different Scenarios

Scenario 1: (Different values on the central lines resistors); here, we consider two different scenarios as in Table 5:
  • 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 6.1) 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 and according to Figure 10 and Table 6, the electricity cost increases for about 30 cents compared to the results obtained in Section 6.1 for the CU profile, since the power loss among the loads and the U G is higher and, as a result, the external grid produces more power to meet the 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 Figure 6b. 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.
Scenario 2: (Faulted lines in the transmission network); here, considering the DU profile, optimal 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.
In the first case, where the transmission line among the loads and the external grid does not work (Figure 11, Table 7), the loss caused by the rest functional lines is higher than in the case of Figure 6b. The UG generated power is distributed through three lines, hence the loss increases as does the electricity cost. A similar situation is observed in the case of a line under fault between the PV and the UG system (Figure 11, Table 7). The power loss is higher towards the UG system from the PV for selling power than in the case of Figure 6b. 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.

7. Conclusions

The goal of this paper was to extend the hierarchical control problem for a meshed DC microgrid already presented in [12] by adding the power losses in the central transmission network. The main objectives were to develop an optimization problem capable of ensuring 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 constrained optimization-based control strategy applying differential flatness, B-spline parametrization and Model Predictive Control methods was introduced.
The mathematical model of the proposed meshed DC microgrid (Figure 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 was obtained in [17]. The DC microgrid was analyzed interpreting all its components, firstly, as electrical circuits. Each component was studied separately, deriving its mathematical model from the associated Bond graph [19,21]. From the Bond graphs, the PH state-space representation [25] of every element was presented in [17]. Next, in this work, the PH dynamical model of the overall central transmission network was provided replacing the transmission lines by resistors.
Then, the hierarchical control approach was analyzed and 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.
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 [26].
The reconfiguration of the system may be accomplished at the high level. The optimization problem developed including the power dissipation within the DC-bus could be used to predict the behavior of the system in the case of unexpected events. 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.
Finally, as mentioned in the paper, the proposed method is applied on a DC microgrid whose parts could be further extended by adding components for energy storage or charges and loads, such as the electrical vehicles, as has been primarily indicated in [17]. The present work gives an analytical study on how to handle the energy management problem and can be further broadened by investigating phenomena like the fast charging, the capacity of electrical transformers to deliver the necessary power and the voltage drops. The proposed method could be also tested on AC networks and components by adding to the existing system different types of AC sources and AC loads.

Author Contributions

Conceptualization: I.Z., I.P., methodology, software, validation, investigation, writing: I.Z.; supervision, I.P., L.L.; All authors have read and agreed to the published version of the manuscript.

Funding

This work is funded by the French National Research Agency within the framework of the project ANR-15-CE05-004-02 C3m (Components, Control and Communication).

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. B-Splines Properties

The following theorem is considered for the B-spline parametrization of the α factor in (23a).
Theorem A1
([27]). Given scalars z ̲ , z ¯ R and a B-spline curve define by z ( t ) = i = 1 n p i b i , d ( t ) and a knot vector T, a sufficient condition for the validation of
z ̲ z ( t ) z ¯
for any t [ τ κ , τ κ + 1 ] T is that:
z ̲ p i z ¯ , κ d + 1 i κ

Appendix B. Supplementary Calculation for the B-Splines

For the constraints (28a), (28b), we give explicitly the calculations for their B-spline parametrization. From (16a), (16b), (25) and the B-splines properties [27], we continue as follows (consider that b i , j , d ( t ) = b i , d ( t ) b j , d ( t ) where 1 i and j n and d = d α ):
v b ( t ) + i b ( t ) R 1 b α ( t ) = i = 1 n p k , i v b b i , d ( t ) + R 1 b i = 1 n p i i b b i , d ( t ) j = 1 n α p j α b j , d α ( t ) = = i = 1 n j = 1 n α ( p i v b + R 1 b p k , i i b ) p j α b i , j , d ( t )
R 1 s c | i b ( t ) | = R 1 s c | i = 1 n p i i b b i , d ( t ) | R 1 s c | i = 1 n p i i b b i , d ( t ) | j = 1 n α p j α b j , d α ( t ) R 1 s c i = 1 n | p i i b | | b i , d ( t ) | j = 1 n α p j α b j , d α ( t ) R 1 s c i = 1 n | p i i b | b i , d ( t ) · · j = 1 n α p j α b j , d α ( t ) = i = 1 n j = 1 n α R 1 s c | p i i b | p j α b i , j , d ( t )
Afterwards, placing (A1a), (A1b) in (29a) and (29b) we obtain:
v D C m i n , h i = 1 n j = 1 n α ( p i v b + R 1 b p i i b ) p j α b i , j , d ( t ) i = 1 n j = 1 n α R 1 s c | p i i b | p j α b i , j , d ( t ) i = 1 n j = 1 n α p i v b + R 1 b p i i b R 1 s c | p i i b | p j α b i , j , d ( t ) ,
v D C m a x , h i = 1 n j = 1 n α ( p i v b + R 1 b p i i b ) p j α b i , j , d ( t ) + i = 1 n j = 1 n α R 1 s c | p i i b | p j α b i , j , d ( t ) i = 1 n j = 1 n α ( p i v b + R 1 b p i i b ) + R 1 s c | p i i b | p j α b i , j , d ( t ) .
where κ d + 1 i , j κ .

References

  1. Roque, L.A.; Fontes, D.B.; Fontes, F.A. A Metaheuristic Approach to the Multi-Objective Unit Commitment Problem Combining Economic and Environmental Criteria. Energies 2017, 10, 2029. [Google Scholar] [CrossRef] [Green Version]
  2. Meng, L.; Sanseverino, E.R.; Luna, A.; Dragicevic, T.; Vasquez, J.C.; Guerrero, J.M. Microgrid supervisory controllers and energy management systems: A literature review. Renew. Sustain. Energy Rev. 2016, 60, 1263–1273. [Google Scholar] [CrossRef]
  3. Ma, J.; Yuan, L.; Zhao, Z.; He, F. Transmission loss optimization-based optimal power flow strategy by hierarchical control for DC microgrids. IEEE Trans. Power Electron. 2016, 32, 1952–1963. [Google Scholar] [CrossRef]
  4. Gamarra, C.; Guerrero, J.M. Computational optimization techniques applied to microgrids planning: A review. Renew. Sustain. Energy Rev. 2015, 48, 413–424. [Google Scholar] [CrossRef] [Green Version]
  5. Iovine, A.; Siad, S.B.; Damm, G.; De Santis, E.; Di Benedetto, M.D. Nonlinear control of a DC microgrid for the integration of photovoltaic panels. IEEE Trans. Autom. Sci. Eng. 2017, 14, 524–535. [Google Scholar] [CrossRef] [Green Version]
  6. Wei, C.; Fadlullah, Z.M.; Kato, N.; Stojmenovic, I. On optimally reducing power loss in micro-grids with power storage devices. IEEE J. Sel. Areas Commun. 2014, 32, 1361–1370. [Google Scholar] [CrossRef]
  7. Feng, K.; Liu, C.; Song, Z. Hour-ahead energy trading management with demand forecasting in microgrid considering power flow constraints. Energies 2019, 12, 3494. [Google Scholar] [CrossRef] [Green Version]
  8. Nahata, P.; La Bella, A.; Scattolini, R.; Ferrari-Trecate, G. Hierarchical Control in Islanded DC Microgrids with Flexible Structures. arXiv 2019, arXiv:1910.05107. [Google Scholar]
  9. Vazquez, N.; Yu, S.S.; Chau, T.K.; Fernando, T.; Iu, H.H.C. A Fully Decentralized Adaptive Droop Optimization Strategy for Power Loss Minimization in Microgrids With PV-BESS. IEEE Trans. Energy Convers. 2018, 34, 385–395. [Google Scholar] [CrossRef]
  10. Sanseverino, E.; Di Silvestre, M.; Badalamenti, R.; Nguyen, N.; Guerrero, J.; Meng, L. Optimal power flow in islanded microgrids using a simple distributed algorithm. Energies 2015, 8, 11493–11514. [Google Scholar] [CrossRef] [Green Version]
  11. Mutarraf, M.; Terriche, Y.; Niazi, K.; Vasquez, J.; Guerrero, J. Energy Storage Systems for Shipboard Microgrids—A Review. Energies 2018, 11, 3492. [Google Scholar] [CrossRef] [Green Version]
  12. Zafeiratou, I.; Prodan, I.; Lefèvre, L.; Piétrac, L. Meshed DC microgrid hierarchical control: A differential flatness approach. Electr. Power Syst. Res. 2020, 180, 106–133. [Google Scholar] [CrossRef]
  13. Zafeiratou, I.; Prodan, I.; Boem, F.; Lefèvre, L. Handling power losses in a DC microgrid through constrained optimization. In Proceedings of the 21st World Congress, Berlin, Germany, 12–17 July 2020. [Google Scholar]
  14. Fliess, M.; Lévine, J.; Martin, P.; Rouchon, P. Differential flatness and defect: An overview. In Workshop on Geometry in Nonlinear Control; Banach Center Publications: Warsaw, Poland, 1993. [Google Scholar]
  15. Lyche, T.; Manni, C.; Speleers, H. Foundations of spline theory: B-splines, spline approximation, and hierarchical refinement. In Splines and PDEs: From Approximation Theory to Numerical Linear Algebra; Springer: Berlin/Heidelberg, Germany, 2018; pp. 1–76. [Google Scholar]
  16. Prodan, I.; Stoican, F.; Louembet, C. Necessary and sufficient LMI conditions for constraints satisfaction within a B-spline framework. In Proceedings of the 2019 IEEE 58th Conference on Decision and Control (CDC), Nice, France, 11–13 December 2019; pp. 8061–8066. [Google Scholar]
  17. Zafeiratou, I.; Prodan, I.; Lefèvre, L.; Piétrac, L. Dynamical modelling of a DC microgrid using a port-Hamiltonian formalism. IFAC-PapersOnLine 2018, 51, 469–474. [Google Scholar] [CrossRef]
  18. Cook, M.D.; Trinklein, E.H.; Parker, G.G.; Robinett, R.D.; Weaver, W.W. Optimal and Decentralized Control Strategies for Inverter-Based AC Microgrids. Energies 2019, 12, 3529. [Google Scholar] [CrossRef] [Green Version]
  19. Duindam, V.; Macchelli, A.; Stramigioli, S.; Bruyninckx, H. Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2009. [Google Scholar]
  20. Paynter, H.M. Analysis and Design of Engineering Systems; MIT Press: Cambridge, MA, USA, 1961. [Google Scholar]
  21. Karnopp, D.C.; Margolis, D.L.; Rosenberg, R.C. System Dynamics: Modeling, Simulation, and Control of Mechatronic Systems; John Wiley & Sons: Hoboken, NJ, USA, 2012. [Google Scholar]
  22. Suryawan, F. Constrained Trajectory Generation and Fault Tolerant Control Based on Differential Flatness and B-Splines. Ph.D. Thesis, School of Electrical Engineering and Computer Science, The University of Newcastle, Callaghan, Australia, 2012. [Google Scholar]
  23. Löfberg, J. YALMIP: A toolbox for modeling and optimization in MATLAB. In Proceedings of the CACSD Conference, Taipei, Taiwan, 2–4 September 2004; Volume 3. [Google Scholar]
  24. Biegler, L.T.; Zavala, V.M. Large-scale nonlinear programming using IPOPT: An integrating framework for enterprise-wide dynamic optimization. Comput. Chem. Eng. 2009, 33, 575–582. [Google Scholar] [CrossRef]
  25. van der Schaft, A.; Jeltsema, D. Port-Hamiltonian systems theory: An introductory overview. Found. Trends Syst. Control 2014, 1, 173–378. [Google Scholar] [CrossRef]
  26. Boem, F.; Ferrari, R.M.; Keliris, C.; Parisini, T.; Polycarpou, M.M. A distributed networked approach for fault detection of large-scale systems. IEEE Trans. Autom. Control 2016, 62, 18–33. [Google Scholar] [CrossRef] [Green Version]
  27. Stoican, F.; Prodan, I.; Popescu, D.; Ichim, L. Constrained trajectory generation for UAV systems using a B-spline parametrization. In Proceedings of the 2017 25th Mediterranean Conference on Control and Automation (MED), IEEE, Valletta, Malta, 3–6 July 2017; pp. 613–618. [Google Scholar]
Figure 1. DC-bus of the DC microgrid with power dissipation included.
Figure 1. DC-bus of the DC microgrid with power dissipation included.
Energies 14 04846 g001
Figure 2. Electrical circuit of ES, where the input voltage v e s is equivalent to the voltage entering from the central transmission network v D C .
Figure 2. Electrical circuit of ES, where the input voltage v e s is equivalent to the voltage entering from the central transmission network v D C .
Energies 14 04846 g002
Figure 3. Bond graph of the DC micogrid’s central transmission network with resistors replacing the transmission lines.
Figure 3. Bond graph of the DC micogrid’s central transmission network with resistors replacing the transmission lines.
Energies 14 04846 g003
Figure 4. Schematic view of the proposed hierarchical control method applied on a meshed DC microgrid.
Figure 4. Schematic view of the proposed hierarchical control method applied on a meshed DC microgrid.
Energies 14 04846 g004
Figure 5. The proposed hierarchical supervision strategy of the DC microgrid considering the power losses in the central transmission network.
Figure 5. The proposed hierarchical supervision strategy of the DC microgrid considering the power losses in the central transmission network.
Energies 14 04846 g005
Figure 6. (a) High level simulation results considering the CU profile. (b) High level simulation results considering the DU profile. Regarding lines in red, they represent the upper and lower part of the corresponding constraints.
Figure 6. (a) High level simulation results considering the CU profile. (b) High level simulation results considering the DU profile. Regarding lines in red, they represent the upper and lower part of the corresponding constraints.
Energies 14 04846 g006
Figure 7. (a) Middle level simulation results considering the CU profile. (b) Middle level simulation results considering the DU profile. Regarding lines in red, they represent the upper and lower part of the corresponding constraints.
Figure 7. (a) Middle level simulation results considering the CU profile. (b) Middle level simulation results considering the DU profile. Regarding lines in red, they represent the upper and lower part of the corresponding constraints.
Energies 14 04846 g007
Figure 8. (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.
Figure 8. (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.
Energies 14 04846 g008
Figure 9. (a) Voltage and current tracking profiles for CU. (b) Voltage and current tracking profiles for DU.
Figure 9. (a) Voltage and current tracking profiles for CU. (b) Voltage and current tracking profiles for DU.
Energies 14 04846 g009
Figure 10. (a) Power balancing and power losses profiles for DU load demand. (b) Power balancing and power losses profiles for CU load demand.
Figure 10. (a) Power balancing and power losses profiles for DU load demand. (b) Power balancing and power losses profiles for CU load demand.
Energies 14 04846 g010
Figure 11. (a) Power balancing and power losses profile generation with R2 equal to 0. (b) Power balancing and power losses profile generation with R4 equal to 0
Figure 11. (a) Power balancing and power losses profile generation with R2 equal to 0. (b) Power balancing and power losses profile generation with R4 equal to 0
Energies 14 04846 g011
Table 1. Parameters.
Table 1. Parameters.
VariableValuesUnits
R 1 s c , R 1 b , 2 b 1, 0.025 , 0.088 [ Ω ]
I 1 s c , I 2 s c 0.25 , 0.25 [H]
C 1 s c , C 2 s c , C 3 s c 0.0008 , 0.0008 , 0.0008 [F]
C 1 b , C 2 b 86 , 400 , 21 , 600 [F]
R 1 , R 2 , R 3 , R 4 1[ Ω ]
Table 2. High and middle level variables and constraints.
Table 2. High and middle level variables and constraints.
VariableValuesUnits
N as in (30a) (30b), (36)27
High level N a as in (30a) (30b), (36)18
d a as in (25)4
v b m i n , h v b v b m a x , h 12 v b 13 [V]
i b m i n , h i b i b m a x , h 9 i b 9 [A]
Constraints q 1 b m i n , m q 1 b q 1 b m a x , m 288.3 q 1 b 307.7 [Ah]
q 2 b m i n , h q 2 b q 2 b m a x , h 72.5 q 2 b 77.5 [Ah]
P u g m i n , h P u g P u g m a x , h 2100 P u g 4200 [W]
v D C m i n , h v u g , p v , e s , l o a d s v D C m a x , h 380 v u g , p v , e s , l o a d s 430 [V]
N p as in (39a)5[h]
Middle level Q v ˜ e s as in (39a) d i a g ( 1 , 1 )
R u ˜ as in (39a)800
v b m i n , m v ˜ b v b m a x , m 11.9 v ˜ b 13.1 [V]
i b m i n , m i ˜ b i b m a x , m 10.6 i ˜ b 10.6 [A]
Constraints q 1 b m i n , m q ˜ 1 b ( k ) q 1 b m a x , m 287.6 q ˜ 1 b ( k ) 308.4 [Ah]
q 2 b m i n , m q ˜ 2 b q 2 b m a x , m 72.3 q ˜ 2 b 77.7 [Ah]
P u g m i n , m P ˜ u g P u g m a x , m 2100 P ˜ u g 4200 [W]
v D C m i n , h v u g , p v , e s , l o a d s v D C m a x , h 370 v u g , p v , e s , l o a d s 430 [V]
Table 3. Generated and consumed power in relation to the corresponding total power.
Table 3. Generated and consumed power in relation to the corresponding total power.
LoadPowerGenerated Power [%]Consumed Power [%]Electricity Cost [Euros]
CU P u g 49.91%1.78% sold to the UG4.318
P e s 0.79%0.76% for battery charging-
P p v 49.30%--
P l o a d s -96.9% for the consumers-
P l o s s -Total: 0.56%
R 1 , 2 , 3 , 4 : 0.12%,
0.24%, 0.11%,
0.09%
-
DU P u g 42.14%13% sold to the UG2.713
P e s 6.58%6.7% for battery charging-
P p v 51.28%--
P l o a d s -79.66% for the consumers-
P l o s s -Total: 0.64%
R 1 , 2 , 3 , 4 : 0.12%,
0.21%, 0.13%,
0.18%
-
Table 4. Generated and consumed power in relation to the corresponding total power.
Table 4. Generated and consumed power in relation to the corresponding total power.
Load ProfilePowerGenerated Power [%]Consumed Power [%]Generated Power Difference from High Level [%]Consumed Power Difference from High Level [%]
CU P u g 49.24%1.85% sold to the UG−0.62%0.07%
P e s 2.08%1.84% for battery charging1.29%1.05%
P p v 48.68%-−0.62%-
P l o a d s -94.74% for the consumers-0.75%
P l o s s -1.57%-1.01%
DU P u g 41.27%13.15% sold to the UG−0.87%0.08%
P e s 7.95%7.73% for battery charging1.37%0.97%
P p v 50.78%-−0.5%-
P l o a d s -77.71% for the consumers-−1.05%
P l o s s -1.41%-0.77%
Table 5. Resistor values for Scenario 1.
Table 5. Resistor values for Scenario 1.
LoadResistorsValueUnit
Commercial R 1 0.5[ Ω ]
R 2 5[ Ω ]
R 3 5[ Ω ]
R 4 0.2[ Ω ]
Domestic R 1 , R 4 1[ Ω ]
R 2 , R 3 3[ Ω ]
Table 6. Simulation results obtained of Scenario 1.
Table 6. Simulation results obtained of Scenario 1.
LoadPowerGenerated Power [%]Consumed Power [%]Electricity Cost (Euros)
CU P u g 51.2%1.9% sold to the UG4.623
P e s 1.9%2% for battery charging-
P l o a d s -92% for the consumers-
P p v 46.9%--
P l o s s -Total: 4.1%
R 1 , 2 , 3 , 4 : 0.19%, 2.1%,
1.76%, 0.05%
-
DU P u g 44%13.1% sold to the UG2.857
P e s 6%6.7% for battery charging-
P l o a d s -78.1% for the consumers-
P p v 50%--
P l o s s -Total: 2.1%
R 1 , 2 , 3 , 4 : 0.2%, 0.9%,
0.6%, 0.4%
-
Table 7. Simulation results obtained of Scenario 2.
Table 7. Simulation results obtained of Scenario 2.
LoadPowerGenerated Power [%]Consumed Power [%]Electricity Cost (Euros)
R 2 = 0 P u g 43%13.1% sold to the UG2.881
P e s 6.8%6.8% for battery charging-
P l o a d s -77.5% for the consumers-
P p v 50.2%--
P l o s s -Total: 2.6%
R 1 , 2 , 3 , 4 : 1.09%, 0%,
0.98%, 0.53%
-
R 4 = 0 P u g 43.7%13.6% sold to the UG2.765
P e s 3.7%3.7% for battery charging-
P l o a d s -80% for the consumers-
P p v 52.6%--
P l o s s -Total: 2.7%
R 1 , 2 , 3 , 4 : 0.98%, 0.61%,
1.11%, 0%
-
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zafeiratou, I.; Prodan, I.; Lefévre, L. A Hierarchical Control Approach for Power Loss Minimization and Optimal Power Flow within a Meshed DC Microgrid. Energies 2021, 14, 4846. https://doi.org/10.3390/en14164846

AMA Style

Zafeiratou I, Prodan I, Lefévre L. A Hierarchical Control Approach for Power Loss Minimization and Optimal Power Flow within a Meshed DC Microgrid. Energies. 2021; 14(16):4846. https://doi.org/10.3390/en14164846

Chicago/Turabian Style

Zafeiratou, Igyso, Ionela Prodan, and Laurent Lefévre. 2021. "A Hierarchical Control Approach for Power Loss Minimization and Optimal Power Flow within a Meshed DC Microgrid" Energies 14, no. 16: 4846. https://doi.org/10.3390/en14164846

APA Style

Zafeiratou, I., Prodan, I., & Lefévre, L. (2021). A Hierarchical Control Approach for Power Loss Minimization and Optimal Power Flow within a Meshed DC Microgrid. Energies, 14(16), 4846. https://doi.org/10.3390/en14164846

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop