A Mathematical Model for the Scheduling of Virtual Microgrids Topology into an Active Distribution Network

: This article presents a method based on a mathematical optimization model for the scheduling operation of a distribution network (DN). The contribution of the proposed method is that it permits the conﬁguration and operation of a DN as a set of virtual microgrids with a high penetration level of distributed generation (DG) and battery energy storage systems (BESS). The topology of such virtual microgrids are modulated in time in response to grid failures, thus minimizing load curtailment, and maximizing local renewable resource and storage utilization as well. The formulation provides the load reduced by bus to balance the system at every hour and the global probability to present energy not supplied (ENS). Furthermore, for every bus, a ﬂexibility load response range is considered to avoid its total load curtailment for small load reductions. The model has been constructed considering a linear version of the AC optimal power ﬂow (OPF) constraints extended for multiple periods, and it has been tested in a modiﬁed version of the IEEE 33-bus radial distribution system considering four different scenarios of 72 h, where the global energy curtailment has been 27.9% without demand-side response (DSR) and 10.4% considering a 30% of ﬂexibility load response. Every scenario execution takes less than a minute, making it appropriate for distribution system operational planning.


Introduction
The concept of resilience has been increasingly used in the literature related to distribution networks (DN) with high penetration levels of distributed energy resources (DER), with an interest in providing greater autonomy and flexibility in the face of unexpected power generation situations [1]. In this context, microgrids (MG) play a relevant role in the transition process to smart grids because they offer an electrical structure that allows for improving the monitoring and management control of the renewable resources on a small scale [2,3]. However, following the definition of MGs given by the U.S Department of Energy [4], a fixed electrical boundary may not be exploiting all the flexibility that an MG provides.
In order to extend the literature, the proposed model addresses this issue and provides additional information for the distribution system operations (DSO). The main contributions of this article are summarized as follows: • A new mathematical optimization model is presented which provides the virtual microgrid topology scheduling of an ADN considering the load curtailment and its flexibility response as the main criteria to partition an unbalanced system into several VMs balanced, with a computational time short enough to ease its adoption by DSO in their usual operational procedures. • We have extended the literature related to the partitioning methods in order to obtain VMs with dynamic boundaries, leaving the algorithms based on hierarchical or graph theory, and providing an optimization linear model. • A quantification of the weaker nodes in terms of the probability of the node to not meet the load under consecutive scenarios, and the size of the energy not supplied which is a useful value in flexibility markets.
The rest of the paper is organized as follows. Section 2 describes the problematic where the mathematical model is situated, and the methodology follows to validate the model under different scenarios. Section 3 presents the mathematical model with the nonlinear constraints already linearized. In Section 4, the model proposed is tested in a case of study and the results are mentioned and discussed. In addition, finally, the conclusions are presented in Section 5.

Preliminaries
This section explains the problem the mathematical model addresses, together with the methodology the authors follow to create a proper test context. Furthermore, a comparison with the current literature is presented.

Problem Definition
Currently, the distribution and transmission systems are supported by different control mechanisms to face the deviation between electricity production and the load [24,25]. This work, typically managed by the transmission system operator (TSO), keeps the energy balancing through three types of control reserves which differ by the temporal space in which they are applied. The question thus arises: why partition a DN to face unexpected events if there are control mechanisms to prevent these kinds of issues? The answer is that the aim of using the resilience concept associated with the DN with high penetration levels of DERs is to dispense with these control mechanisms and be able to achieve self-adaptation to a failure of high impact on the grid [26]. Therefore, the ADN must be able to measure the energy amount that cannot be met due to the absence of the utility grid and to identify the load points which could be applied to a load curtailment in order to minimize the energy losses.
The DSR concept is related to the amount of energy which the ADN could reduce to avoid an unbalanced system under the absence of the utility grid [27]. i.e., if the ADN possesses the flexibility to reduce a percentage of its energy load, the global curtailment will be lower than the case of an ADN where there is no chance to reduce the load more than to zero to avoid an unbalanced system. However, the flexibility degree of the DSR must be known in advance by the ADN such that the model can decide when and how much load to reduce.

Related Works
The literature related to the optimal partition of an ADN, introduced above, is grouped into algorithm or optimization models to explain the main differences between the current methodologies and the formulation proposed.
The work presented by [14] makes an overview of the clustering techniques available in the literature, together with a brief summary of its applicability. Some of these algorithms have been used to partition a DN into MGs or VMs, for which the most relevant are those who use community detection [8,[17][18][19][20] based on a different index, such as: electrical distance [19,20], electrical coupling strength [17], or voltage sensitivity matrix [18] because they allow an identification, through an index, of the weakest point of a DN and the use of these weakest points as limit points between clusters. These methodologies are powerful tools when the partition required is for fixed boundaries and the technical network constraints are considered. However, the limitation appears when the topology is variable through time and some new elements, such as BESS and DSR, are introduced into the ADN since these flexibility elements possess a variable impact that must be considered in a time horizon where the next step depends on the behavior of the previous step.
The optimization models are an alternative of the algorithms because they permit finding an optimal value and include elements external to the DN, such as the BESS and DSR. Nevertheless, their main drawbacks lay in the nonlinearity and non-convexity of the AC-OPF formulation, which significantly increases their computation burden [28,29] and prevent their direct application for multi-period problems in an operational context. Specifically, the MINLP includes, in its AC version, the cosine and sine functions, which are the main problem because the current solvers available in the market possess high execution times. Under this context, some linearization methods are applied to reduce the complexity of the AC-OPF [30][31][32][33], of which some have been used in clustering optimization problems. However, there are few publications that use optimization models for the scheduling of VMs. The work in [23] uses a linearization and bender decomposition method to tackle the formation of dynamic MGs for an ADN with a high penetration of DERs. However, it does this without considering a flexible response of the loads. On the other hand, [22] uses spectral graph theory in a bi-level optimization problem to find MG reconfiguration into a DN with 123-buses, but the formulation does not include storage systems and DSR, and it is not applied in a multi-period context. Therefore, considering the sparse literature on the topic, one of the goals of the model proposed here is to extend the literature related to the optimization models applied in a clustering method to form VMs, considering DSR to increase the resilience of an ADN, and regardless of the current external mechanism.

Methodology
The methodology followed in this article consists of two stages: the first one assumes the capacity and location that should have the DG and BESS into the DN. The second stage addresses the capability of the DN to face unexpected failures through the self-partition into VMs, which is the main contribution of this article.
The placement and sizing have been addressed separately in [13], where a novel algorithm presented allows finding an optimal solution to the capacity and location for different types of DGs and BESSs, considering a time horizon longer than 24 h and basing the formulation on an AC-OPF to include the technical constraints of the network (this is represented in the upper level of Figure 1). The lower level of the figure corresponds to the optimal topology of the DN in case the power energy is insufficient to meet the electricity load. Therefore, following Figure 1, the optimal model for the VMs topology scheduling (VMTS) is executed for scenarios of 72 h, where every hour has two possible results regarding the system: the first one being, where the capacity installed, weather conditions, and the electricity load are balanced, meaning that the DN does not need to be partitioned into VMs. However, for the second result, if the electricity load is not satisfied by the power generation of the DGs, the model finds the amount of energy load that must be reduced at each bus to balance the system, where every load point possesses a range of load that can be reduced. If the load curtailment exceeds this range, the load in that point is reduced to zero and the bus is isolated, coming in a bound of VM. Thus, for each hour where these conditions appear, the model provides the scheduling of the topology that must follow the DN to continuously balance the system, sacrificing, of course, a load curtailment that is minimizing on the objective function but avoiding an outage.
In order to make sure that the topologies found by the mathematical model are feasible for the classic AC-OPF problem, the loads associated with the isolated nodes are replaced by zeros and the buses are modified with a reduced load within its range of flexibility. The DN with the new set of loads is executed in a multi-period AC-OPF. Thus, the feasible solutions of the VMTS model are also feasible for the classic AC-OPF problem.

Optimal Model for Scheduling Virtual Microgrids Topology
The formulation explained in this article is a modified version of the model presented in our previous work [13], which uses as the main structure an AC-OPF extended to a multi-period and a linearization done by [28]. However, unlike [13], which finds the optimal place and location of DERs into a DN, the model here is focused on finding the scheduling of a dynamic topology of an ADN to face lower generation scenarios or failure on the utility grid, through the minimization of the load curtailment in order to find the weakest nodes, and identify eventual flexibilities. The model is explained in two subsections; the first one corresponds to the base structure to find the proper topology of the ADN to face unbalance issues without taking into account a flexible demand response, and a second section, where the power consumption is not a parameter anymore but a variable, in order to include flexibility on demand.

Model without Demand Flexibility
The virtual microgrids topology scheduling (VMTS) model is described as follows: The objective function in Equation (1) minimizes the active (δ p i,t ) and reactive (δ q i,t ) power not supplied, and the active (p t i,j ) and reactive (q t i,j ) power flow between the nodes. Thus, this structure allows us to: • Promote the fulfillment of the electricity load since the positive variables δ p i,t and δ q i,t are penalized (Υ) when their values are greater than zero, if and only if, the power generation and the energy storage is not enough to meet the load. • Prioritize the loads with or closer to a DG or BESS when it is not a chance to meet the power consumption through the minimization of the flow variables (p t i,j ) and (q t i,j ).
In order to penalize the flow between nodes, the absolute value for the active and reactive power flow between nodes is considered, which can take negative and positive values depending on the flow direction, and because an absolute value requires a simple change of variable to become linear. On the other hand, the ENS penalization (Υ) is a parameter that is defined by the DSO regarding its power quality, i.e., if the ENS is more relevant than the energy losses through the lines, then Υ is greater than the power flow through the lines. In the opposite case, Υ must be equal to or lower than the power flow in lines.
The formulation constraints are as follows: Equations (2) and (3) correspond to active and reactive nodal power flow equations; therefore, ∀i ∈ B, ∀t ∈ T, Both equations, and the remaining constraints, have already been linearized using the methodology presented in [28], holding the left side of these equations as the classic AC-OPF. However, on the right side, we consider the active power injected by the distributed resource (pg k i,t ); a virtual active power (λ p i,t ), which represents the power remaining to meet the load in case of low generation; the active power demanded by the node (PL i,t ), wherein this formulation is a parameter; and the power absorbed (ch l i,t ) or injected (ds l i,t ) from the battery to the grid. The same structure follows Equation (3).
The set of constraints from (4) to (7) correspond to the linear version of the active and reactive power flow equation which have been widely explained in previous AC-OPF formulations [13]; therefore, ∀(i, j) ∈ L, ∀t ∈ T, Power load and virtual power are explained in constraints (8)-(11); therefore, ∀i ∈ B, ∀t ∈ T, Equation (8) represents the total power load which depends on the power supplied (γ p i,t ) and not supplied (δ p i,t ) in every node i, thus, the total load is composed by these two types of variable, where the power not supplied is penalized in the objective function, to ensure that δ p i,t is as small as possible. Equation (9), the power not-delivered (δ p i,t ), bound the virtual power λ p i,t which represents the remaining power that it is needed to meet the total power load (see Equation (2). Therefore, if the power generation is insufficient to provide the power required, the virtual power (λ p i,t ) variable fulfills this gap, being penalized in the objective function, if and only if this gap is equal to or lower than the power not supplied (δ p i,t ) in the node i. Constraints (10) and (11) describe the same logic except for the reactive power load.
Active and reactive power limits together with the voltage security constraints are considered in Equations (12)-(14); therefore, ∀i ∈ B, ∀k ∈ G type , Equation (12) makes sure that the voltage in each node and at every period is within the limits allowed. The parameters are squared due to the change of variable needed to linearize the original power flow equation. Constraints (13) establish the limits for the power injected by the DG in the node through two different parameters; U g l corresponds to the total power capacity installed and PG max k,t is the power generated by the DG in a certain period, which is a normalized power generation curve. Therefore, Equation (13) establishes the maximum power that can be injected in Equation (2) at every hour, where PG max k,t follows the power generation behavior of a photovoltaic panel (PV)/wind turbine (WT), and it allows for including the stochastic behavior of the DG into the OPF constraints.
The sets of constraints (15)-(18) represent a linear version of the apparent power limit; therefore, where: Equations (15) and (16) emerge as a set of constraints which represent a polygon formed by r ∈ R lines that define a feasible region approximation of the original apparent power limit, which is a circle. Thus, the bigger the set R, the closer the polygon will be to the circle.
To model the behavior of the batteries, the following constraints are added. Therefore, ∀i ∈ B, ∀t ∈ T, ∀l ∈ B type , The state of charge (SOC) for battery type l is defined in Equation (19) and is composed of the SOC of the previous period plus the energy injected and absorbed from the system with their respective efficiency factors. In order to maintain the structure of the formulation of our previous work [13], the SOC is in MWh and not in percentage, as usual; thus, the SOC represented in Equation (19) is bounded by SOC max l /SOC min l which are in % and are multiplied by the total capacity installed. In this way, the soc l i,t belongs only to the capacity range operation of the battery. Therefore, for a capacity of 10 MWh and a 90%/10% as operational limits, the soc l i,t+1 will not be less than 1 MWh or greater than 9 MWh. On the other hand, the maximum power injected or taken from the DN for a certain period is explained in constraints (21) and (22) through M l which is multiplied by a binary variable representing when the battery is charging or discharging. In the case that the battery is discharging, i.e., w l i,t = 0, the maximum energy discharged must be equal to or lower than M l because the parameter Y l i is one when there is a battery in the node. Otherwise, if Y l i (a binary parameter) is zero, the constraints (23) are zero and all the Equations (19)-(22).

Model with Demand Flexibility
In this section, Equations (1), (2), (8) are rewritten in order to include the DSR into the model explained previously, and a new constraint is added to bound the flexibility degree of the energy demand: Equation (24) is the new objective function for which, unlike (1), the amount of power consumption reduced (µ i,t ) is minimized together with the energy flow and the active/reactive power not supplied, in order to balance the system. However, this new variable µ i,t must be lower than or equal to the demand response of every load. Therefore, ∀i ∈ B, ∀t ∈ T where Ω i,t in Equation (25) represents the percentage of demand response of every load at every time, i.e., this vector parameter represents the flexibility of every load to reduce its demand. On this way, different flexibilities can be considered. For example, Ω 1,5 = 0 means that the bus 1 at 5th hour cannot reduce its demand. On the other hand, Ω 2,5 = 0.15 means that bus 2 at the 5th hour can reduce its demand by 15%. In addition, µ i,t must be integrated into Equations (2) and (8) to keep the balance in the power flow equations, where the parameter PL i,t passes to be an initial power consumption where the flexibility response is subtracted. Thus, ∀i ∈ B, ∀t ∈ T

Results and Discussion
In this section, the two versions of the mathematical model are applied under four different scenarios of 72 hours. Moreover, the impact of the DSR on the VMs formed by the model proposed is presented and discussed.

Case Study
A multi-period version of the PG&E 33-bus radial distribution network [34] has been used to test the model proposed. The sizing and location of the DERs have been separately obtained using the algorithm presented in our previous work [13], considering a connected-grid system because an autonomous network, i.e., the energy, is provided 100% by DERs; it is not a proper case of study to observe the capability of the ADN to form sub-systems when the grid fails. Therefore, using [13], the location and sizing have been set as shown in Figure 2.
In order to test the mathematical optimization model, four different 72-hour scenarios corresponding to each season of the year are taken into account to study the different topologies formed by the ADN when the system is disconnected from the utility grid and the weather conditions are not capable of supplying 100% of the power consumption. The weather resources are taken from historical data (a year and a time-step of 1[h]) of the Lugo, Galicia region in the north of Spain 2018; Ref.
[35] for the wind speed, and [36] for the temperature and irradiance. The power generation data have been obtained using the modeling presented in [37] for the wind turbines, and the formulation used in [38] for PV panels. Thus, the energy demand and the power generation for the periods of study are presented in Figure 3. Every scenario of 72 hours has been chosen from the 8760 h available in the data set, with the aim to represent periods of 72 h with a high difference between both the power consumption and generation. Moreover, every curve has been normalized to respect the order of magnitude and comment on the power capability available. For example, the fourth period/scenario possesses a high-level wind resource, which means that the wind turbine will be able to produce energy at its maximum capacity. On the other hand, there are many consecutive hours that show a very low power capability for PV panels and wind turbines, which will force load curtailment and the partitioning of the ADN into sub-systems. The DSR parameter is set at 30% for every load at every time, i.e., every bus at every hour has the chance to reduce its load by 30%.
The mathematical model has been programmed in Python using Cplex as the main solver, and it has been run on an Intel Core i7-10510U CPU @ 1.80GHz 2.30 GHz with 16 GB memory and operative system Windows, 64 bit.  It should be noted that the goal of the case study is; (1) to size the DSR effect on the DN topology scheduling, and (2) to show the type of data provided by the model and how this information can be useful in terms of operational planning for DSO. It is important to mention this factor because the results depend exclusively on the DERs capacity, electricity demand, season, geographic location, weather conditions, and the DSR considered.

Model Performance
The mathematical formulation, explained in Section 3, has been compared with its MINLP version to observe the huge time differences between nonlinear and linear versions. Thus, Table 1 shows the times of execution for both models and their respective objective functions (OF). When the model is not multi-period, both formulations reach the same OF. However, when the model is extended to multi-periods, the computation burden for the MINLP increases faster than the MILP to the point that at hour 12 the time for the MINLP is around 5.5 h, which is much longer than the alternative version that takes only 3.79 s, keeping an error close to the 7%. Due to the exponential time execution increase for the MINLP version, Table 1 has been completed until hour 12 and not beyond. Nevertheless, the MILP running time for 72 h is close to a minute that puts it within the operational range of DSO.

Results
In order to improve the understanding of the data obtained from the model, the results are grouped as follows: • Model operation explains how the model meets the energy load through the available resources and shows the role of the storage devices.
• The weakest points identify the nodes with a high probability of ENS. This probability considers any percentage of non-compliance from the total load for a certain bus at every hour. • Flexibility Market presents the amount of ENS for every node at every hour and compares the total ENS when DSR is considered and when it is not. • Dynamic topologies show the topology scheduling that must follow the ADN to pass from an unbalanced system to a balanced system, with lower load compliance than the original.

Mathematical Model Operation
Figures 4 and 5 are presented to show how the model proposed to operate to meet the energy demand with the available resources ( Figure 4) and observe the impact of the storage devices ( Figure 5) for a certain scenario without considering DSR. Thus, the fourth period has been considered where, in the first hours, it is possible to meet the energy load with the DERs, due to existing high wind resources and irradiance. Furthermore, the SOC of the batteries increases during the first half of the scenario, achieving its maximum, to face the second half where the weather conditions are not the best. From the 40th hour on, the batteries are discharged (see Figure 5) and the first nodes with low energy not supplied by the system start to appear because the discharge rate of the battery does not allow for fulfilling the gap between load and generation. Finally, at the end of the period, there is not enough energy stored in the batteries to meet the energy demand. To improve the visualization of Figure 4, the energy provided by the DERs has been overlapped to observe the share of the energy demand that is satisfied by each resource, and, for Figure 5, the curves have been normalized, taking as reference the maximum value of each curve.
In global terms, for the 288 h considered, composed of four scenarios of 72 h, the total energy demanded by the system is 434.17 MWh/pr, of which 56.3% is provided by WT, 8.2% by PV panels, 13.4% the storage systems, and 4.6% of energy losses. The ENS when the DSR is not taken into account is 27.9%; however, when 30% of DSR is considered, i.e., every load point can reduce its energy load at a maximum of 30%; otherwise, the load is reduced to zero, the ENS drop-down 62.7%, from 27.9% to 10.4%. Therefore, for the installed capacity and for the considered scenarios, the system could provide 73.3% of the energy as maximum, and the remaining energy load would be reduced or traded through an eventual energy flexibility market.

The Weakest Nodes
The weak nodes of the ADN have been represented in Table 2 through the probability of every bus not meeting the total electricity load in a period of 72h. i.e., for a scenario of 72 h, if, in every hour of this scenario, a certain bus does not supply the entire energy demand, and its probability of default is 100%. The first column indicates the buses, where bus 1 has not been considered because it does not have a load at any time. The following four columns indicate the different scenarios of 72 h considered, and the column Total shows the probability of every bus for ENS considering a 288 h period, i.e., assuming the four scenarios as a set of consecutive hours. The remaining columns possess the same meaning, but considering the DSR version model, and the last column DERs installed indicates whether the bus possesses or is next to a DER. Therefore, following the results presented in the total columns, the buses with a DG or BESS installed have a lower probability of default when the weather conditions are not proper to generate enough power to meet demand with the capacity installed, and the support of the utility grid is not guaranteed. Specifically, the buses with WT installed are the ones that have the lowest probability of default followed by the nodes with PV panels and the nodes with a storage system. Moreover, as the objective function (1) and (24) minimize the energy flow through the lines when there is non-compliance, the nodes which are furthermost of the buses with a DER system, are the points with the highest probability to present load curtailment, as is the case of the upper section of the DN test system {23, 24, 25} and the lower section {19, 20, 21, 22}.
The first two scenarios present the highest probability of ENS with or without DSR. The results are expected because they are the scenarios with more consecutive hours with low power generation (see Figure 3), in contrast with scenarios 3 and 4 which have a better power generation available. However, even though the ENS trend is the same, the probability decreases significantly in some buses when the DSR is included in the model. The clearest example is the case of scenario 4, where the ENS probability is zero for every bus when the DSR is considered.

Energy for Flexibility Market
Unlike Table 2, which shows the probability of default, Table 3 presents the amount of energy that it is not possible to fully supply to the system, for every bus, and for every scenario. For example, in the fourth scenario, for the second bus in Table 3, it would be required to generate an extra 0.11 MWh to meet the electricity load on the weather conditions considered, and not produce a load curtailment for the second bus, in order to avoid an unbalanced system. Therefore, this extra energy for each node and for every hour could be used by the DSO to be traded in an eventual energy flexibility market. On the other hand, for the same previous example, when the DSR is considered into the formulation, the extra 0.11 MWh can be reduced (because it is within the flexibility range response) and therefore a total load shaved is not needed. In addition, the right part of Table 3 shows a clear decrement of ENS due to the DSR effect, which is visualized easily in the fourth scenario, which does not present load reduced. The proportion of ENS for every node follows the same pattern of Table 2 because, for the nodes with power generation or storage devices, the amount of extra energy needed to fulfill the electricity load is lower than the buses without DERs.

Dynamic Topology
The dynamic topology offered by the model is through the variable associated with the ENS (δ p i,t ), which takes values greater than zero when the DERs cannot provide sufficient energy to meet the demand. These values are recorded for every period at each bus. Table 4 shows simplified results of these values, where the values greater than the maximum DSR flexibility have been replaced by one, and when these values are within the range allowed by the DSR or the load is satisfied, the values are zero for the sake of clarity. Specifically, Table 4 represents the results for the time comprised between period 40 h to 48 h corresponding to the second scenario. This scenario has been selected because it presents many buses with load curtailment, which is the case at hour 46. Taking as a reference the results in Table 4, hour 46 has been depicted in Figure 6 to observe how the VMs, represented by binary numbers in the previous table, take form in the ADN. The red points indicate the nodes where the load has been decreased to zero, and the dashed red lines to visualize the isolated system which is bounded by the red points. For those nodes where the load has been suppressed but there is power generation, the energy injected in these red points is used to meet the demand in the closest points of the ADN for which the load is not zero. On the other hand, when the DSR is included in the formulation, the load curtailment is reduced, as shown on the right side of Table 4, leaving only buses 21 and 22 without energy supplied and formed the VM represented in Figure 7.

Discussion
The application of the results presented above focuses mainly on the operation planning of a DN. The mathematical model provides the following five key points to face unexpected failures in the utility grid: • Simulate a failure scenario of 72 h. • Low computational burden due to the linearity of the formulation. • Obtain the virtual microgrids topology scheduling that should follow the DN to keep the system balanced at all times, considering a DSR factor. • Quantify the amount of energy that must be reduced or traded an eventual energy flexibility market, at every bus, and at every hour. • Identify the probability of every bus to present a load curtailment, which sizes the service level under the failure scenario.
Nevertheless, there are some limitations that could emerge if the model is run for scenarios greater than 168 h or if the model is applied for DN with a high density of buses because the time response of the model could exceed the time response of the DSO for the operation planning. Therefore, future research could focus on testing the computational burden of the model for different test distribution systems greater than 33-buses.

Conclusions
This paper has presented an optimization model to find the topology that must be followed for an ADN to self-balance under unexpected grid failures, identifying the load reduction required to balance the system again, and isolate the loads, with a curtailment greater than its capacity response, which become the limits of the VMs formed. Every new topology has been validated using the classic AC-OPF model, in order to make sure that the feasible solution of the model proposed is also a feasible solution for the model used by the DSO.
The model has been tested under four scenarios of 72 h to consider the seasonality of renewable resources, and the location and sizing of the DERs have been taken into account as parameters known. Therefore, the nodes which have presented a high probability of ENS and a high load curtailment have been further away from the DG and BESS. This result is explained because the objective function minimizes the flow between the nodes so that the model prioritizes the loads closer to a DG or BESS. In addition, when DSR is included in the load parameters, the global load curtailment is reduced by 62.7%, i.e., from 27.9% when the DSR is not included in the mathematical model, to a 10.4% when the loads are capable of self-modified in order to avoid an unbalanced system. New research related to this field can take place, introducing the constraints presented in this article within the location and sizing problems for DG and BEES into DN, to make sure that the placement obtained by the model takes into account the load curtailment in the case of grid failures.