Mitigating Energy System Vulnerability by Implementing a Microgrid with a Distributed Management Algorithm

: This work presents a management strategy for microgrid (MG) operation. Photovoltaic (PV) and wind generators, as well as storage systems and conventional units, are distributed over a wide geographical area, forming a distributed energy system, which is coordinated to face any contingency of the utility company by means of its isolated operation. The management strategy divides the system into three main layers: renewable generation, storage devices, and conventional units. Interactions between devices of the same layer are determined by solving an economic dispatch problem (EDP) in a distributed manner using a consensus algorithm (CA), and interactions between layers are determined by means of a load following strategy. In this way, the complex behaviour of PV and wind generation, the battery storage system, and conventional units has been effectively combined with CA to solve EDP in a distributed manner. MG performance and its vulnerability are deeply analysed by means of an illustrative case study. From the observed results, vulnerability under extreme conditions could be reduced up to approximately 30% by coupling distributed renewable generation and storage capacity with an energy system based on conventional generation.


Introduction
The constant development of information and communication technologies and the depletion of natural resources accompanied by climate change have exposed infrastructures required for industrial and manufacturing processes and the provision of services and products to new, unique threats. As a consequence, the definition of vulnerability and the manners to mitigate it are being widely studied.
Energy systems are undergoing an important transformation by adopting communication capabilities to improve monitoring procedures and consumer interactions. This represents a transition from a traditional energy system to an intelligent one.
Similar to other industrial infrastructures, an intelligent system is exposed to risks from the interaction between the power system and the communication infrastructure (CI), as well as natural phenomena [1]. A representative example is the December 2015 cyber-attack in Ukraine, where a six-hour blackout affected hundreds of thousands of consumers in Kiev [2].
To provide a resilient power supply, this panorama has encouraged the development of the microgrid (MG) concept. In this regard, problems caused by Hurricane Sandy on the East Coast of the United States have motivated the installation of an MG as a part of the energy system rebuilding plan [3,4]. Similarly, an MG installed at Stafford Hill is able to provide energy to critical infrastructures in a resilient and effective manner [5].
Resilient power generation can be achieved by linking together equipment spread over a community to build an MG capable of operating in stand-alone mode. Thus, distributed energy Load profile of small commercial loads [19].
The methodology used for load modelling is based on [20] and has been widely used for the analysis of HESs. It aims to create a yearly time series from the daily profile of Figure 2. The first step consists of creating a correlated random time series according to (1). Randomness is introduced by using the variable ( ) , which is modeled as a Gaussian random number with a mean equal to zero and standard deviation equal to √1 − (Ø ) 2 . The second step consists of normalising the daily profile of Figure 2. The variable ( , ) is obtained by repeating the daily profile of interest ( Figure 2) in a periodic manner. Then, the normalised time series ( ( , ) ) is obtained by applying (2). The addition of correlated ( ( , ) ) and normalised ( ( , ) ) time series according to (3) results in a time series with some correlation degree (Ø ) and the diurnal profile of Figure 2. Based on the equality presented in (4), which is in fact a probability transformation, the simulated load time series ( ( , ) ) can be obtained according to (5). The cumulative distribution function (CDF) of the simulated load demand ( ) can be modeled by using a log-normal distribution or beta distribution [21].

Model of Distributed PV Generator
A mathematical model of the distributed PV generator is presented in (7)- (23). The whole system is composed of several ( ) PV generators. For a determined PV system ( ), cell temperature ( ( , ) ) is shown in (7) [20]. Equations (8)- (16) [22] represent the model of the PV cell, which establishes the relationship between voltage, current, and power. Simulated maximum power point tracking (MPPT) has been performed by applying the golden section search algorithm (GSSA) on (16) according to [23]. The rated capacity of the power converter has been estimated using (17), while the Load profile of small commercial loads [19].
The methodology used for load modelling is based on [20] and has been widely used for the analysis of HESs. It aims to create a yearly time series from the daily profile of Figure 2. The first step consists of creating a correlated random time series according to (1). Randomness is introduced by using the variable ξ (t) , which is modeled as a Gaussian random number with a mean equal to zero and standard deviation equal to 1 − (Ø L ) 2 . The second step consists of normalising the daily profile of Figure 2. The variable L (n,t) is obtained by repeating the daily profile of interest ( Figure 2) in a periodic manner. Then, the normalised time series (l (n,t) ) is obtained by applying (2). The addition of correlated (ϕ (n,t) ) and normalised (l (n,t) ) time series according to (3) results in a time series with some correlation degree (Ø L ) and the diurnal profile of Figure 2. Based on the equality presented in (4), which is in fact a probability transformation, the simulated load time series (P T(n,t) ) can be obtained according to (5). The cumulative distribution function (CDF) of the simulated load demand ( f L ) can be modeled by using a log-normal distribution or beta distribution [21]. ϕ (n,t) = Ø L ϕ (n,t−1) + ξ (t) ; (1) l (n,t) = L (n,t) − µ σ ; (2) z (n,t) = l (n,t) + ϕ (n,t) ; (3) f Z z (n,t) = f L P T(n,t) ; (4) P L(n,t) =   P T(n,t) − min P T(n,t) max P T(n,t) − min P T(n,t)   P max L − P min L + P min L .
Characteristics of demand time series (P T(n,t) ) can be modified by scaling the values between two different limits (P min L and P max L ) according to (6), resulting in a time series (P L(n,t) ) with a diurnal profile similar to that shown in Figure 2, limited within a pre-defined load interval.

Model of Distributed PV Generator
A mathematical model of the distributed PV generator is presented in (7)- (23). The whole system is composed of several (P) PV generators. For a determined PV system (p), cell temperature (T PV(p,t) ) η I(p,t) = p I(p,t) p I(p,t) + α I(p) + β I(p) p I(p,t) + θ I(p) p I(p,t) 2 .

Model of Distributed Wind Generator
Power production from a typical wind generator can be estimated by means of a power curve similar to the one shown in Figure 3 [20]. Equations (24)-(26) present the cost curve for each wind generator. In constraint (27), P max W(w,t) is directly obtained from the evaluation of Figure 3, and P opt W(w,t) is obtained from the implementation of CA.

Model of Distributed Wind Generator
Power production from a typical wind generator can be estimated by means of a power curve similar to the one shown in Figure 3 [20]. Equations (24)-(26) present the cost curve for each wind generator. In constraint (27), ( , ) is directly obtained from the evaluation of Figure 3, and ( , ) is obtained from the implementation of CA.

Model of Distributed BESS
Due to its promising characteristics, a vanadium redox battery (VRB) was selected among BESS technologies. The mathematical model previously proposed in [25][26][27] was used in its modelling. The simulation approach consists of scaling a single battery of 5 kW/20 kWh by connecting elements in serial (N S B(b) ) and parallel (N P B(b) ). Important variables, such as voltage, power, efficiency, and state of charge (SOC) during charging (P B(b,t) > 0) and discharging (P B(b,t) ≤ 0) conditions have been estimated. Definitions of battery voltage, efficiency, and SOC are shown in (28), (29), and (30), respectively, while some of their operating constraints are presented in (31) and (32). Calculation of battery voltage during charging and discharging is explained in (33) and (34), respectively. The general values for voltage and energy efficiency for charging processes are computed according to (35)-(37); for discharging processes, these variables can be calculated by using (38)-(40). Constraint (41) indicates that battery power has to be limited to the rating power of the cell stack. The operating cost of BESS is described in (42)-(44), which represents the cost of cycling energy through the storage system. Constraint (45) is related to the power dispatch for charging and discharging processes (P opt B(b,t) ); therefore, the power to be charged or discharged to or from BESS has to be limited to the charging and discharging power available (P B(b,t) ) at a determined time instant t. Values of the parameters involved on VRB modelling are reported in Table 1.

Model of Distributed CG
CGs are capable of providing controllable power within a determined interval. Power generation is adjusted to minimise the operating cost, which is described according to (46) and rewritten in (47) and (48). Power generation should be limited to its minimum (P min C(g) ) and maximum (P max C(g) ) operating values according to (49).

Architecture and CI
As stated previously, the management strategy proposed in this paper for vulnerability mitigation on DSs is based on EDP solved by CA. Figure 4 describes the general structure of CI, where all PV and wind generators are connected to the same infrastructure. Similarly, BESS and CGs interact within their corresponding communications network. The energy management system (EMS) continuously interacts with the DS through the smart grid so that distributed HES can start operation under any contingency. During normal operating conditions, HES is supposed to be fully integrated on the electricity market to maximise its profits by managing renewable generators, BESS, and CGs in the appropriate manner. The energy management system (EMS) continuously interacts with the DS through the smart grid so that distributed HES can start operation under any contingency. During normal operating conditions, HES is supposed to be fully integrated on the electricity market to maximise its profits by managing renewable generators, BESS, and CGs in the appropriate manner.
Once the distributed autonomous operation has been activated to face any threat, the EMS distributes PV and wind generation with the highest priority to take advantage of renewable sources. Distributed BESS provides storage capabilities for excess renewable resources, which can be used later when no wind or solar energy is available. Finally, fossil-fuel-based generation is used to supply the power needs not yet satisfied. CA implementation proposed by Yang et al. [28] has been used in this paper to carry out the solution of EDP for distributed PV and wind generation, BESS, and CG. CA is traditionally formulated to solve EDP in systems vertically integrated, considering the cost curves previously described in (20)  Once the distributed autonomous operation has been activated to face any threat, the EMS distributes PV and wind generation with the highest priority to take advantage of renewable sources. Distributed BESS provides storage capabilities for excess renewable resources, which can be used later when no wind or solar energy is available. Finally, fossil-fuel-based generation is used to supply the power needs not yet satisfied.
CA implementation proposed by Yang et al. [28] has been used in this paper to carry out the solution of EDP for distributed PV and wind generation, BESS, and CG. CA is traditionally formulated to solve EDP in systems vertically integrated, considering the cost curves previously described in Instead of considering the operating cost for each technology, the energy selling price specified by the generation/storage device owner has been adopted. In the case of PV generation, δ 1 PV(p) → 0 , δ 2 PV(p) > 0, and δ 3 PV(p) = 0 have been assumed, where δ 2 PV(p) is the energy selling price of the corresponding PV generator. Similarly, δ 1 W(w) → 0 , δ 2 W(w) > 0, and δ 3 W(w) = 0 have been assumed for each wind generator. δ 28 B(b) → 0 , δ 29 B(b) > 0, and δ 30 B(b) = 0 have been assumed for each BESS, while δ 1 C(g) → 0 , δ 2 C(g) > 0, and δ 3 C(g) = 0 have been assumed for each CG. This approach incorporates the economic benefit of each generator and storage device.
As mentioned, PV and wind generators are integrated into a single CI managed by means of CA. For the sake of simplicity, CA implementation only considering PV generation will be described. For the other generation/storage technologies, the formulation is similar. Minimum (λ min PV(p) ) and maximum (λ max PV(p) ) incremental cost values are shown in (50) and (51), respectively. Evolution of incremental cost ( Energies 2019, 12, 616 10 of 30 . . .
Equations (50)-(57) can also be used to estimate the power dispatch of systems based on wind generation, BESS, and CG. Assuming P min C(g) → 0 , CG owners can offer power production at a determined selling price (δ 2 C(g) ) limited to its rated capacity (P max C(g) ). Under this context, the proprietary of each device of the energy system (PV and wind generator, BESS, and CG) only needs to specify its energy price and power available.

MG Management Strategy
As this work aims to create a management strategy to reduce DS vulnerability, the primary effort has been applied to analysis of the system under emergency. Variables involved in the system modelling and management are shown in Figure 5. EMS is continuously monitoring state variables of the DS, such as voltage (U S(n,t) ) and power (P S(n,t) and Q S(n,t) ) at medium voltage (MV). The operation of EMS is performed at a low voltage (LV) level using local resources. The management technique is based on a load following strategy [15], wherein renewable generation and BESS are committed to reducing fuel-based power sources and increasing the long-term autonomy of energy systems.
Specifically, two operating modes have been considered. In the first, MG is operated to prevent voltage problems. In the second, MG operates independently from the utility company, thus avoiding total collapse. These modes are carefully described in the next subsections.
Energies 2018, 11, x FOR PEER REVIEW 10 of 29 operation of EMS is performed at a low voltage (LV) level using local resources. The management technique is based on a load following strategy [15], wherein renewable generation and BESS are committed to reducing fuel-based power sources and increasing the long-term autonomy of energy systems. Specifically, two operating modes have been considered. In the first, MG is operated to prevent voltage problems. In the second, MG operates independently from the utility company, thus avoiding total collapse. These modes are carefully described in the next subsections.

Microgrid Operation to Avoid Voltage Collapse
Depending on the technical design of the DS under study, the rate of demand growth, lack of investment in adequacy and overhauling of the DS, and other factors, problems related to voltage stability could occur during extreme operating conditions. In this regard, Venkatesh et al. [29] proposed the evaluation of MLI as a quantitative measure of how far the system is from voltage collapse. MLI depends on active and reactive power flows at the MV feeder, as well as the system voltage. MLI is evaluated according to (58): Voltage problems occur when the condition ( , ) < 1 is fulfilled. This condition is then resolved by reducing the power flow on ∆ ( , ) , which is calculated using (59):

.Microgrid Operation under Substation Contingency
Distributed PV and wind generators are tightly related to distributed-BESS operation. As previously explained, MG management is based on a load following strategy, wherein the surplus of renewable generation is used for charging VRB. However, the amount of power being purchased from renewable generation must be known to operate the entire MG cost-effectively. In other words, the amount of power not consumed by the MG energy demand and not stored on VRB can be understood as a cost increment; however, this condition can be prevented by carefully predetermining the power to be supplied by renewable sources. Figure 6 presents the algorithm

Microgrid Operation to Avoid Voltage Collapse
Depending on the technical design of the DS under study, the rate of demand growth, lack of investment in adequacy and overhauling of the DS, and other factors, problems related to voltage stability could occur during extreme operating conditions. In this regard, Venkatesh et al. [29] proposed the evaluation of MLI as a quantitative measure of how far the system is from voltage collapse. MLI depends on active and reactive power flows at the MV feeder, as well as the system voltage. MLI is evaluated according to (58): Voltage problems occur when the condition MLI (n,t) < 1 is fulfilled. This condition is then resolved by reducing the power flow on ∆P MLI(n,t) , which is calculated using (59): ∆P MLI(n,t) = MLI (n,t) − 1 P S(n,t) + jQ S(n,t) . (59)

Microgrid Operation under Substation Contingency
Distributed PV and wind generators are tightly related to distributed-BESS operation. As previously explained, MG management is based on a load following strategy, wherein the surplus of renewable generation is used for charging VRB. However, the amount of power being purchased from renewable generation must be known to operate the entire MG cost-effectively. In other words, the amount of power not consumed by the MG energy demand and not stored on VRB can be understood as a cost increment; however, this condition can be prevented by carefully predetermining the power to be supplied by renewable sources. Figure 6 presents the algorithm applied to dispatch the available PV and wind power. If the available renewable power is higher than the MG load, the power surplus should be used for charging VRB. Under these conditions, the VRB maximum charging power is calculated by using (28)-(33), (35)-(37), and (41). This procedure results in generating the power required by the entire BESS (∑ P B ), which is later added to the MG load (P L ) to estimate the power to be supplied by PV and wind systems. If the available renewable power is higher than the power required for charging VRB and for satisfying the MG demand, CA is applied according to Section 2.6 to calculate the power to be purchased from each PV (P opt PV(p,t) ) and wind system (P opt W(w,t) ), taking into account their respective energy prices (δ 2 PV(p) and δ 2 W(w) ). On the contrary, if renewable resources are not enough to supply the MG demand and VRB charging power, they are entirely committed without any optimization process. In other words, P opt PV(p,t) ← P PV(p,t) and P opt W(w,t) ← P W(w,t) are assumed.
Under these conditions, the VRB maximum charging power is calculated by using (28)-(33), (35)-(37), and (41). If excess energy results in lower than the maximum energy to be stored by BESS, EDP by CA is applied to determine how energy charging is distributed. The energy distribution among BESS depends upon their respective prices. Conversely, if the renewable power dispatched is not enough to supply the energy demand, BESS is fully discharged. The VRB maximum discharging power is calculated by using (28)-(32), (34), and (38)-(41). CG provides the remaining power to supply the MG demand. The flowchart presented in Figure  8 describes the management algorithm. If distributed CG is capable of supplying the entire demand, power dispatch is performed by using CA. In the contrary case, all units are dispatched at their maximum or rated capacity, and the load demand is cut to some degree. Once the contribution of renewable sources (P opt PV(p,t) and P opt W(w,t) ) has been determined, the results are used to analyse the behaviour of distributed VRB. Management of BESS is carried out using the flowchart shown in Figure 7. Depending on the amount of renewable energy available, BESS is charged or discharged.  If dispatched renewable power is higher than the load demand, distributed VRB is charged. Under these conditions, the VRB maximum charging power is calculated by using (28)-(33), (35)-(37), and (41). If excess energy results in lower than the maximum energy to be stored by BESS, EDP by CA is applied to determine how energy charging is distributed. The energy distribution among BESS depends upon their respective prices. Conversely, if the renewable power dispatched is not enough to supply the energy demand, BESS is fully discharged. The VRB maximum discharging power is calculated by using (28)-(32), (34), and (38)-(41).
CG provides the remaining power to supply the MG demand. The flowchart presented in Figure 8 describes the management algorithm. If distributed CG is capable of supplying the entire demand, power dispatch is performed by using CA. In the contrary case, all units are dispatched at their maximum or rated capacity, and the load demand is cut to some degree.

Case Studies
The proposed approach to reduce the vulnerability of a smart DS is illustrated in this section by analysing a hypothetical energy system to be located in Puerto Rico (latitude 18.022º and longitude -66.0168º). Time series of wind speed, solar radiation ( ( ) ), and ambient temperature ( ( ) ) were found in [18]. A load demand time series was created according to the methodology described in Section 2 using the parameters specified in Table 2 and a time step of 1 h (∆ = 1 h). Data of distributed PV and wind generators are presented in Tables 3 and 4, respectively. The number of PV cells on each panel ( ( ) ) and the ideality factor ( ) were assumed to be 1. BESS and CG data are shown in Tables 5 and 6, respectively. The proposed management strategy was implemented in MATLAB © on a personal computer with an i7-3630QM CPU at 2.4 GHz with 8 GB of memory and a 64-bit operating system.

Case Studies
The proposed approach to reduce the vulnerability of a smart DS is illustrated in this section by analysing a hypothetical energy system to be located in Puerto Rico (latitude 18.022 • and longitude −66.0168 • ). Time series of wind speed, solar radiation (G (t) ), and ambient temperature (T A(t) ) were found in [18]. A load demand time series was created according to the methodology described in Section 2 using the parameters specified in Table 2 and a time step of 1 h (∆t = 1 h). Data of distributed PV and wind generators are presented in Tables 3 and 4, respectively. The number of PV cells on each panel (N CP PV(p) ) and the ideality factor (m PV ) were assumed to be 1. BESS and CG data are shown in Tables 5 and 6, respectively. The proposed management strategy was implemented in MATLAB© on a personal computer with an i7-3630QM CPU at 2.4 GHz with 8 GB of memory and a 64-bit operating system.    The nominal voltage of the DS was assumed to be 4.16 kV. The peak load was assumed to be 700 kVA, with a power factor equal to 95%. The distribution feeder was assumed to be built using a conductor 26.24 kcmil (7 strands) type AAC, operating in a system of 60 Hz. CA was considered as having a tolerance value of 0.001 (ε = 0.001) and 150 iterations (I = 150). CI for the management of renewable power sources, BESS, and CG are shown in Figures 9-11, respectively. These structures were used to build the weighting matrices for each layer of Figure 4 (CI PV P and CI PV Q for PV generators). The following subsections describe in detail how the proposed management strategy can be used to avoid voltage collapse on a typical distribution feeder and how an MG can be operated when the substation is unable to provide service.   The following subsections describe in detail how the proposed management strategy can be used to avoid voltage collapse on a typical distribution feeder and how an MG can be operated when the substation is unable to provide service.

Case Study 1: MG Operation to Avoid Voltage Collapse
In this subsection, the operation of distributed HES to prevent voltage collapse of a distribution feeder is analysed. For the sake of simplicity, the MG is composed of only CGs with the characteristics previously shown in Table 6 and CI presented in Figure 11. As explained in subsection 3.1, MLI was evaluated using (58). Then, at those hours when ( , ) < 1, power from CGs was dispatched to mitigate the instability of the distribution feeder. As mentioned, the power of CGs (∆ ( , ) ) is adjusted according to (59). Figure 12 presents load demand time series taking into account the effects of CGs. System instability was observed 128 times, which represents 1.461187% of the year duration (loss of load probability equal to 1.461187%).   The following subsections describe in detail how the proposed management strategy can be used to avoid voltage collapse on a typical distribution feeder and how an MG can be operated when the substation is unable to provide service.

Case Study 1: MG Operation to Avoid Voltage Collapse
In this subsection, the operation of distributed HES to prevent voltage collapse of a distribution feeder is analysed. For the sake of simplicity, the MG is composed of only CGs with the characteristics previously shown in Table 6 and CI presented in Figure 11. As explained in subsection 3.1, MLI was evaluated using (58). Then, at those hours when ( , ) < 1, power from CGs was dispatched to mitigate the instability of the distribution feeder. As mentioned, the power of CGs (∆ ( , ) ) is adjusted according to (59). Figure 12 presents load demand time series taking into account the effects of CGs. System instability was observed 128 times, which represents 1.461187% of the year duration (loss of load probability equal to 1.461187%).   The following subsections describe in detail how the proposed management strategy can be used to avoid voltage collapse on a typical distribution feeder and how an MG can be operated when the substation is unable to provide service.

Case Study 1: MG Operation to Avoid Voltage Collapse
In this subsection, the operation of distributed HES to prevent voltage collapse of a distribution feeder is analysed. For the sake of simplicity, the MG is composed of only CGs with the characteristics previously shown in Table 6 and CI presented in Figure 11. As explained in subsection 3.1, MLI was evaluated using (58). Then, at those hours when ( , ) < 1, power from CGs was dispatched to mitigate the instability of the distribution feeder. As mentioned, the power of CGs (∆ ( , ) ) is adjusted according to (59). Figure 12 presents load demand time series taking into account the effects of CGs. System instability was observed 128 times, which represents 1.461187% of the year duration (loss of load probability equal to 1.461187%).

Case Study 1: MG Operation to Avoid Voltage Collapse
In this subsection, the operation of distributed HES to prevent voltage collapse of a distribution feeder is analysed. For the sake of simplicity, the MG is composed of only CGs with the characteristics previously shown in Table 6 and CI presented in Figure 11. As explained in Section 3.1, MLI was evaluated using (58). Then, at those hours when MLI (n,t) < 1, power from CGs was dispatched to mitigate the instability of the distribution feeder. As mentioned, the power of CGs ∆P MLI(n,t) is adjusted according to (59). Figure 12 presents load demand time series taking into account the effects of CGs. System instability was observed 128 times, which represents 1.461187% of the year duration (loss of load probability equal to 1.461187%). For a better understanding, Figure 13 shows the first week (168h) of time series previously shown in Figure 12. As can be observed, the MG operation reduces the load demand during peak load hours. In other cases, this problem is mitigated by means of demand response resources. In this paper, the operation of properly coordinated distributed resources has been used.  Figure 14 shows the power production of each CG for the first week; such a power dispatch is optimised according to the prices specified by the owner of each generator (Table 6).  For a better understanding, Figure 13 shows the first week (168h) of time series previously shown in Figure 12. As can be observed, the MG operation reduces the load demand during peak load hours. In other cases, this problem is mitigated by means of demand response resources. In this paper, the operation of properly coordinated distributed resources has been used. For a better understanding, Figure 13 shows the first week (168h) of time series previously shown in Figure 12. As can be observed, the MG operation reduces the load demand during peak load hours. In other cases, this problem is mitigated by means of demand response resources. In this paper, the operation of properly coordinated distributed resources has been used.  Figure 14 shows the power production of each CG for the first week; such a power dispatch is optimised according to the prices specified by the owner of each generator (Table 6).   Figure 14 shows the power production of each CG for the first week; such a power dispatch is optimised according to the prices specified by the owner of each generator (Table 6).
Probability distribution of MLI for both series, previously shown in Figure 12 (with and without MG operation), is presented in Figure 15. According to these results, the incorporation of an MG system based on distributed CG guarantees the stable functioning of a DS. When CGs are committed to supplying power to maintain system stability (∆P MLI(n,t) ), MLI takes a value of 1 or higher. Probability distribution of MLI for both series, previously shown in Figure 12 (with and without MG operation), is presented in Figure 15. According to these results, the incorporation of an MG system based on distributed CG guarantees the stable functioning of a DS. When CGs are committed to supplying power to maintain system stability (∆ ( , ) ), MLI takes a value of 1 or higher. For example, Table 7 and Figures 16-18 show the results for = 37 h. Power dispatch is shown in Table 7, where it can be observed how required power can be cost-effectively provided. Evolution of CA is shown in Figure 16, which converges to the optimal values shown in Table 7. Power mismatch and incremental cost are shown in Figures 17 and 18, respectively.   Probability distribution of MLI for both series, previously shown in Figure 12 (with and without MG operation), is presented in Figure 15. According to these results, the incorporation of an MG system based on distributed CG guarantees the stable functioning of a DS. When CGs are committed to supplying power to maintain system stability (∆ ( , ) ), MLI takes a value of 1 or higher. For example, Table 7 and Figures 16-18 show the results for = 37 h. Power dispatch is shown in Table 7, where it can be observed how required power can be cost-effectively provided. Evolution of CA is shown in Figure 16, which converges to the optimal values shown in Table 7. Power mismatch and incremental cost are shown in Figures 17 and 18, respectively.  For example, Table 7 and Figures 16-18 show the results for t = 37 h. Power dispatch is shown in Table 7, where it can be observed how required power can be cost-effectively provided. Evolution of CA is shown in Figure 16, which converges to the optimal values shown in Table 7. Power mismatch and incremental cost are shown in Figures 17 and 18, respectively.

Case Study 2: Long-Term Autonomous MG Operation
In this subsection, the operation of distributed HES to prevent a power shortage of critical loads is analysed. In this case, the MG operates independently of the DS. Simulations were performed with a time step of 1 h (∆t = 1 h). The load-demand time series used in this case is shown in Figure 19.

Case Study 2: Long-Term Autonomous MG Operation
In this subsection, the operation of distributed HES to prevent a power shortage of critical loads is analysed. In this case, the MG operates independently of the DS. Simulations were performed with a time step of 1 hour (∆ = 1 h). The load-demand time series used in this case is shown in Figure 19.   Figure 21 presents VRB power, and Figure 22 shows VRB SOC. On the right side of Figure 21, it is possible to observe how the cheapest BESS ( = 5) has the most important role due to its low selling price and its considerable high storing capacity ( (5) = 20). It is also possible to observe how the second BESS ( = 2) gets the highest SOC due to its low storing capacity ( (2) =10) and its relatively low selling price.

Case Study 2: Long-Term Autonomous MG Operation
In this subsection, the operation of distributed HES to prevent a power shortage of critical loads is analysed. In this case, the MG operates independently of the DS. Simulations were performed with a time step of 1 hour (∆ = 1 h). The load-demand time series used in this case is shown in Figure 19.   Figure 21 presents VRB power, and Figure 22 shows VRB SOC. On the right side of Figure 21, it is possible to observe how the cheapest BESS ( = 5) has the most important role due to its low selling price and its considerable high storing capacity ( (5) = 20). It is also possible to observe how the second BESS ( = 2) gets the highest SOC due to its low storing capacity ( (2) =10) and its relatively low selling price.  Figure 21 presents VRB power, and Figure 22 shows VRB SOC. On the right side of Figure 21, it is possible to observe how the cheapest BESS (b = 5) has the most important role due to its low selling price and its considerable high storing capacity (N P B(5) = 20). It is also possible to observe how the second BESS (b = 2) gets the highest SOC due to its low storing capacity (N P B(2) =10) and its relatively low selling price.   Finally, Figure 23 shows time series of a CG power dispatch, which illustrates how power generation is similarly shared around the cheapest units.

Case Study 3: Short-Term Autonomous MG Operation
In this subsection, short-term operation of the distributed HES is studied. In this case, a time step of 15 minutes (∆ = 0.25 h) was considered. Power dispatch on an hourly basis requires a forecasting tool in four steps ahead of 15 minutes each. A forecasting tool can be based on time series analysis [30] so that historical data recorded by EMS can be used to accomplish new tasks.  Finally, Figure 23 shows time series of a CG power dispatch, which illustrates how power generation is similarly shared around the cheapest units.

Case Study 3: Short-Term Autonomous MG Operation
In this subsection, short-term operation of the distributed HES is studied. In this case, a time step of 15 minutes (∆ = 0.25 h) was considered. Power dispatch on an hourly basis requires a forecasting tool in four steps ahead of 15 minutes each. A forecasting tool can be based on time series analysis [30] so that historical data recorded by EMS can be used to accomplish new tasks. Finally, Figure 23 shows time series of a CG power dispatch, which illustrates how power generation is similarly shared around the cheapest units.  Finally, Figure 23 shows time series of a CG power dispatch, which illustrates how power generation is similarly shared around the cheapest units.

Case Study 3: Short-Term Autonomous MG Operation
In this subsection, short-term operation of the distributed HES is studied. In this case, a time step of 15 minutes (∆ = 0.25 h) was considered. Power dispatch on an hourly basis requires a forecasting tool in four steps ahead of 15 minutes each. A forecasting tool can be based on time series analysis [30] so that historical data recorded by EMS can be used to accomplish new tasks.

Case Study 3: Short-Term Autonomous MG Operation
In this subsection, short-term operation of the distributed HES is studied. In this case, a time step of 15 min (∆t = 0.25 h) was considered. Power dispatch on an hourly basis requires a forecasting tool in four steps ahead of 15 min each. A forecasting tool can be based on time series analysis [30] so that historical data recorded by EMS can be used to accomplish new tasks.
For illustrative purposes, information publically available at [31] has been scaled according to the magnitudes previously presented in Tables 3 and 4 for a typical day. In other words, renewable power and load time series of 24 h with a step of 0.25 h were artificially created for simulation purposes; each series has 96 values. Figure 24 shows load time series of MG under analysis. It is supposed to be the result of a determined forecasting process applied to get a one-hour prediction (four steps ahead) that works 24 h of a typical day. The load series was created by scaling to a maximum value of 700 kVA and power factor close to 1. Figure 25 presents PV and wind generation time series artificially created as previously mentioned. For illustrative purposes, information publically available at [31] has been scaled according to the magnitudes previously presented in Tables 3 and 4 for a typical day. In other words, renewable power and load time series of 24 h with a step of 0.25 h were artificially created for simulation purposes; each series has 96 values. Figure 24 shows load time series of MG under analysis. It is supposed to be the result of a determined forecasting process applied to get a one-hour prediction (four steps ahead) that works 24 hours of a typical day. The load series was created by scaling to a maximum value of 700 kVA and power factor close to 1. Figure 25 presents PV and wind generation time series artificially created as previously mentioned.   Figures 26 and 27 show the simulation of VRB in terms of power and SOC, respectively. According to these results, the VRB operation is highly influenced by the behaviour of renewable generation, especially PV.
Initially, VRB should be discharged to reduce the load demand jointly with wind generation. Excess PV generation leads VRB to be continuously charged. It can be observed that the relationship between the energy prices and storing capacities is similar to those previously discussed in subsection 4.2 (Case study 2).  For illustrative purposes, information publically available at [31] has been scaled according to the magnitudes previously presented in Tables 3 and 4 for a typical day. In other words, renewable power and load time series of 24 h with a step of 0.25 h were artificially created for simulation purposes; each series has 96 values. Figure 24 shows load time series of MG under analysis. It is supposed to be the result of a determined forecasting process applied to get a one-hour prediction (four steps ahead) that works 24 hours of a typical day. The load series was created by scaling to a maximum value of 700 kVA and power factor close to 1. Figure 25 presents PV and wind generation time series artificially created as previously mentioned.   Figures 26 and 27 show the simulation of VRB in terms of power and SOC, respectively. According to these results, the VRB operation is highly influenced by the behaviour of renewable generation, especially PV.
Initially, VRB should be discharged to reduce the load demand jointly with wind generation. Excess PV generation leads VRB to be continuously charged. It can be observed that the relationship between the energy prices and storing capacities is similar to those previously discussed in subsection 4.2 (Case study 2).   Figures 26 and 27 show the simulation of VRB in terms of power and SOC, respectively. According to these results, the VRB operation is highly influenced by the behaviour of renewable generation, especially PV.  Finally, Figure 28 shows how CGs work at similar power production due to their similarities in electricity prices (Table 6) established by their corresponding owners.  Finally, Figure 28 shows how CGs work at similar power production due to their similarities in electricity prices (Table 6) established by their corresponding owners. Initially, VRB should be discharged to reduce the load demand jointly with wind generation. Excess PV generation leads VRB to be continuously charged. It can be observed that the relationship between the energy prices and storing capacities is similar to those previously discussed in Section 4.2 (Case study 2).
Finally, Figure 28 shows how CGs work at similar power production due to their similarities in electricity prices (Table 6) established by their corresponding owners.

Vulnerability Assessment of Autonomous MG
In this subsection, the vulnerability [32] of MG configuration of Figure 4 is analysed. Specifically, CIs shown in Figures 9 and 11 have been used to evaluate the vulnerability of the whole system against contingencies intentionally created (CIC) on the renewable and CG subsystems, respectively. Failures or contingencies in the distributed BESS layer have not been considered in this study.
In this regard, Figures 29 and 30 explain how the contingencies were treated. First, vertex 1 of each CI was assumed to be the command vertex, which receives information from EMS. Through this vertex, the load demand at each time instant is reported to each layer.
Once any edge is removed from CI, those vertices directly connected to the command vertex (vertex 1) are assumed to be active and consequently under the control of EMS. Contrarily, those vertices out of the compass of the command vertex are supposed to be offline (inactive). After all CICs are realized, information is exchanged among active vertices forming a ring. This is done because a strongly connected directed graph is required for solving EDP by CA. Finally, EDP is solved by following the procedure explained in Section 3, just considering the active vertices.
At the left side of Figure 29, edges between vertices 2 and 4, 4 and 7, 6 and 7, and 5 and 6 have been removed from the CI of renewable generation. Under this scenario, vertices 4, 5, and 6 are inactive, while vertices 1, 2, 3, 7, 8, and 9 remain connected to the command vertex (vertex 1). At the right side of Figure 29, it is observed how the information is changed according to a ring topology, as previously explained.

Vulnerability Assessment of Autonomous MG
In this subsection, the vulnerability [32] of MG configuration of Figure 4 is analysed. Specifically, CIs shown in Figures 9 and 11 have been used to evaluate the vulnerability of the whole system against contingencies intentionally created (CIC) on the renewable and CG subsystems, respectively. Failures or contingencies in the distributed BESS layer have not been considered in this study.
In this regard, Figures 29 and 30 explain how the contingencies were treated. First, vertex 1 of each CI was assumed to be the command vertex, which receives information from EMS. Through this vertex, the load demand at each time instant is reported to each layer.

Vulnerability Assessment of Autonomous MG
In this subsection, the vulnerability [32] of MG configuration of Figure 4 is analysed. Specifically, CIs shown in Figures 9 and 11 have been used to evaluate the vulnerability of the whole system against contingencies intentionally created (CIC) on the renewable and CG subsystems, respectively. Failures or contingencies in the distributed BESS layer have not been considered in this study.
In this regard, Figures 29 and 30 explain how the contingencies were treated. First, vertex 1 of each CI was assumed to be the command vertex, which receives information from EMS. Through this vertex, the load demand at each time instant is reported to each layer.
Once any edge is removed from CI, those vertices directly connected to the command vertex (vertex 1) are assumed to be active and consequently under the control of EMS. Contrarily, those vertices out of the compass of the command vertex are supposed to be offline (inactive). After all CICs are realized, information is exchanged among active vertices forming a ring. This is done because a strongly connected directed graph is required for solving EDP by CA. Finally, EDP is solved by following the procedure explained in Section 3, just considering the active vertices.
At the left side of Figure 29, edges between vertices 2 and 4, 4 and 7, 6 and 7, and 5 and 6 have been removed from the CI of renewable generation. Under this scenario, vertices 4, 5, and 6 are inactive, while vertices 1, 2, 3, 7, 8, and 9 remain connected to the command vertex (vertex 1). At the right side of Figure 29, it is observed how the information is changed according to a ring topology, as previously explained.  Vulnerability assessment is performed by means of hourly simulations. In this regard, yearly time series of wind speed, solar radiation, ambient temperature, and load demand are used. The study consists of estimating the load lost as a result of several CIC. A CIC consists of removing a determined edge or edges randomly chosen from a determined CI.
The procedure for vulnerability assessment can be explained with greater details through the following steps: First, obtain the hourly (∆ = 1 h) time series of renewable resources and load demand; second, select a determined number of edges to be removed from a determined CI; third, vertices out of the compass of the command vertex are supposed to be offline (inactive). After all CICs are realized, information is exchanged among active vertices forming a ring. This is done because a strongly connected directed graph is required for solving EDP by CA. Finally, EDP is solved by following the procedure explained in Section 3, just considering the active vertices.
At the left side of Figure 29, edges between vertices 2 and 4, 4 and 7, 6 and 7, and 5 and 6 have been removed from the CI of renewable generation. Under this scenario, vertices 4, 5, and 6 are inactive, while vertices 1, 2, 3, 7, 8, and 9 remain connected to the command vertex (vertex 1). At the right side of Figure 29, it is observed how the information is changed according to a ring topology, as previously explained.
At the left side of Figure 30, edges between vertices 2 and 3, 4 and 5, 5 and 6, and 6 and 7 have been removed from the CI of the CG. Under this context, vertices 3, 4, 6, and 7 are inactive, while vertices 1, 2, and 5 remain connected to the command vertex (vertex 1). At the right side of Figure 30, it is observed how the information in changed according to a ring topology.
Vulnerability assessment is performed by means of hourly simulations. In this regard, yearly time series of wind speed, solar radiation, ambient temperature, and load demand are used. The study consists of estimating the load lost as a result of several CIC. A CIC consists of removing a determined edge or edges randomly chosen from a determined CI.
The procedure for vulnerability assessment can be explained with greater details through the following steps: First, obtain the hourly (∆t = 1 h) time series of renewable resources and load demand; second, select a determined number of edges to be removed from a determined CI; third, considering those vertices connected to the command vertex, as well as the renewable resources and load demand at each time instant (t), solve EDP by CA; fourth, determine the load lost (kW) for each hour of the year; fifth, calculate the load lost in percentage by dividing the annual load lost and the MG annual load; and sixth, calculate the connectivity loss of CI by dividing the number of edges to be removed (Step 2) and the total amount of edges on the CI under study. This procedure is repeated until all edges on CI are removed.
Considering the topology of Figure 29, if the number of edges to be removed is set to 5, it means that five edges randomly chosen would be removed at each hour (t = 1, . . . , T = 8760 h) of the year. Then, the loss of load would be estimated in percentage from the yearly simulation, and connectivity loss would be equal to 41.6% (5/12). This procedure is repeated until 100% connectivity loss is reached.
The vulnerability of an energy system used in Section 4.1 (distributed CG system only) was compared to the vulnerability of that used in Section 4.2 (distributed PV/wind/BESS/CG system) when only contingencies on a CG subsystem were considered. This exercise revealed the importance of renewable generation and energy storage on the vulnerability of a conventional system. Corresponding results are presented in Figure 31, where it can be observed how renewable generation and energy storage reduce the load to be supplied by a CG, which reduces the load lost and, consequently, the vulnerability of the whole system. On the other hand, under extreme conditions (connectivity loss equal to 100%, as in Figure 31), the vulnerability reduces by approximately 30%.
Similarly, vulnerability analysis considering contingencies on renewable generation and CG subsystems is presented in Figure 32. This analysis incorporates the results previously reported in Figure 31 with additional information. When connectivity loss is equal to 0% on a renewable generation subsystem, the vulnerability curve shown in Figure 31 for the distributed PV/wind/BESS/CG system (blue line) is obtained. From this point, it is possible to observe how fast the whole system collapses when the renewable subsystem is subject to intentional contingencies. This fact reveals the important interdependence between renewables and CG subsystems. Similarly, vulnerability analysis considering contingencies on renewable generation and CG subsystems is presented in Figure 32. This analysis incorporates the results previously reported in Figure 31 with additional information. When connectivity loss is equal to 0% on a renewable generation subsystem, the vulnerability curve shown in Figure 31 for the distributed PV/wind/BESS/CG system (blue line) is obtained. From this point, it is possible to observe how fast the whole system collapses when the renewable subsystem is subject to intentional contingencies. This fact reveals the important interdependence between renewables and CG subsystems.

Conclusions and Remarks
In this paper, HES to reduce DS vulnerability has been studied, and detailed computational models representing behaviours of renewable generation, VRB, and fossil fuel power generation have been combined with a decentralised management strategy to develop an energy system able to face the negative consequence of utility company malfunctioning. Management of the proposed configuration has been discussed under two different conditions: voltage stability of the DS and autonomous operation of the entire distributed HES.
Traditionally, HESs located in remote areas are frequently composed of single power generation and storage devices. In addition, these devices belong to a single owner so that they are managed to  Similarly, vulnerability analysis considering contingencies on renewable generation and CG subsystems is presented in Figure 32. This analysis incorporates the results previously reported in Figure 31 with additional information. When connectivity loss is equal to 0% on a renewable generation subsystem, the vulnerability curve shown in Figure 31 for the distributed PV/wind/BESS/CG system (blue line) is obtained. From this point, it is possible to observe how fast the whole system collapses when the renewable subsystem is subject to intentional contingencies. This fact reveals the important interdependence between renewables and CG subsystems.

Conclusions and Remarks
In this paper, HES to reduce DS vulnerability has been studied, and detailed computational models representing behaviours of renewable generation, VRB, and fossil fuel power generation have been combined with a decentralised management strategy to develop an energy system able to face the negative consequence of utility company malfunctioning. Management of the proposed configuration has been discussed under two different conditions: voltage stability of the DS and autonomous operation of the entire distributed HES.
Traditionally, HESs located in remote areas are frequently composed of single power generation and storage devices. In addition, these devices belong to a single owner so that they are managed to

Conclusions and Remarks
In this paper, HES to reduce DS vulnerability has been studied, and detailed computational models representing behaviours of renewable generation, VRB, and fossil fuel power generation have been combined with a decentralised management strategy to develop an energy system able to face the negative consequence of utility company malfunctioning. Management of the proposed configuration has been discussed under two different conditions: voltage stability of the DS and autonomous operation of the entire distributed HES.
Traditionally, HESs located in remote areas are frequently composed of single power generation and storage devices. In addition, these devices belong to a single owner so that they are managed to minimise the generating cost, satisfying load requirements. In the configuration considered in this paper, elements of HES are dispersed over a wide region composed of many power generation and storage devices. Each of these elements can be owned by a different investor, which specifies its appropriate energy selling price. Moreover, HES for power supply in isolated areas is designed to provide energy in a cost-effective manner, while the distributed configuration proposed in this paper is designed to provide electricity for minimising system vulnerability by splitting power generation and storage devices into many components. Of note, investments in DS updating and overhauling can be directed to increase the deployment of distributed systems around the geographical regions with load demand of high relevance, not just concentrating economic resources on DS. The results observed from the analysis of the case studies illustrate how the proposed management strategy can be used to perform the coordination of generators and BESSs as independent aggregated units, where the owner of each device establishes a determined price by the energy produced (generation units) or cycled (storage devices). Moreover, vulnerability analysis revealed how important interdependence is between renewable and conventional generation. In fact, a reduction of 30% for system vulnerability was observed with respect to the traditional system only based on CG.

Conflicts of Interest:
The authors declare no conflict of interest.

T PV(p,t)
Photovoltaic cell temperature of generator p at time t ( • C). G (t) Solar radiation at time t (W/m 2 ). NOCT (p) Nominal operating cell temperature of generator p ( • C).  σ Standard deviation of load demand time series (kW). Ø L Autocorrelation coefficient.
ε Parameter of consensus algorithm implementation. l (n,t) Normalized load demand at time t for node n (kW). L (n,t) Load demand at time t for node n (kW). z (n,t) Correlated and profiled load demand at time t for node n. f Z (·) Cumulative distribution of z (n,t) time series.
f L (·) Cumulative distribution of P T(n,t) time series. → f PV (·) Projection operator for dispatch of photovoltaic system.

P S(n,t)
Flow of active power of branch n at time t (kW). P T(n,t) Simulated load demand at time t for node n (kW). P L(n,t) Scaled load demand at time t for node n (kW). P min L Minimum load demand (kW). P max L Maximum load demand (kW). P PV(p,t) Photovoltaic cell power at time t for generator p (W). P I(p) Rated power of converter of generator p (kW). p I(p) Normalized power of inverter at time t for generator p.    Optimal power dispatch of battery system b at time t (kW).

PV(p)
Parameters of cost curve of photovoltaic generator p.

W(w)
Parameters of cost curve of wind generator w.

C(g)
Parameters of cost curve of conventional generator g.
η PV(p) Cell efficiency of generator p.
η I(p,t) Inverter efficiency at time t for generator p. Battery discharging power efficiency at time h. F PV(p,t) Cost curve of conventional generator p at time t ($). F W(w,t) Cost curve of conventional generator w at time t ($). Cost curve of conventional generator b at time t ($). F C(g,t) Cost curve of conventional generator g at time t ($). Resistance of branch n (Ω).

X S(n)
Reactance of branch n (Ω). MLI (n,t) Maximum loadability index at time t for node n.