Review on Energy Resilience

Contents lists available at ScienceDirect

Applied Energy

journal homepage: www.elsevier.com/locate/apenergy

A linear two-stage method for resiliency analysis in distribution systems considering renewable energy and demand response resources

Saeed Mousavizadeh⁎, Mahmoud-Reza Haghifam, Mohammad-Hossein Shariatkhah Faculty of Electrical and Computer Engineering, Tarbiat Modares University, Tehran, PO Box 14115-111, Iran

H I G H L I G H T S

• The concept of the resiliency in distribution networks is addressed. • A linear path-based approach is proposed to model the topological characteristics of the network. • A two-stage stochastic programming approach is applied to characterize the uncertainty of the renewables. • Impacts of the demand side programs and renewable energies on system resiliency are investigated.

A R T I C L E I N F O

Keywords: Resiliency Distribution network Microgrids Linear programming Distributed energy resource

A B S T R A C T

Due to the unique structure and special characteristics of electric distribution networks, along with the increase in the number and severity of the natural disasters in recent years, presenting a proper framework and procedure for evaluating the resiliency of these systems is more essential. In addition, it is obligatory to model the different capabilities and features of the smart grids components in order to provide a better knowledge about their expected performance in such circumstances. In response to these challenges, this paper addresses the concept of resiliency and its dimensions in distribution networks. A new model based on mixed-integer linear programming is proposed to properly model and evaluate the resiliency of smart distribution systems. In the proposed model, optimal formation of dynamic microgrids (MGs), their service areas, and the optimal management of different technologies such as energy storage (ES) units, demand side management programs and distributed generations (DGs) units are investigated. In addition, employing a two-stage framework based on stochastic programming, the impact of increasing penetration level of the renewable energy resources and their related uncertainties on system resiliency is examined. The efficiency and applicability of the proposed integrated model is verified by performing multiple simulations on the modified 118-bus test system and a real distribution network.

1. Introduction

Due to the climatic changes, the number and intensity of the natural disasters such as hurricanes, floods, droughts, etc. have increased sig- nificantly in many countries in recent years [1–3]. Enlargement of the destructive effects of such phenomena on electric power system implies that the utilization of the existing methods and common strategies for the operation and planning of power systems which have been devel- oped based on the familiar concepts such as the risk, reliability, stabi- lity, security, etc. will be faced with serious challenges in the near fu- ture [4,5]. In such incidents, even systems with high reliability and appropriate status in terms of risk and equipment vulnerability may encounter significant problems and experience extensive blackouts. Hence, in recent years, the concept of resiliency as a complement to the

previous concepts is introduced to ensure the optimal functionality of the power system in these circumstances [6,7]. A resilient system will be able to efficiently implement the necessary actions to reduce the effects of the catastrophic events [8].

Up until now, the issue of power system resiliency has been ex- amined in several papers, and some indices have been developed for this purpose [9–11]. However, owing to the higher vulnerability of the distribution network against these events, and its special structure and distinct characteristics, providing suitable frameworks to evaluate and enhance the system behavior in such situations will be of great im- portance. Despite the large body of research in this field, less attention has been paid to resiliency concept and its dimension in distribution networks.

Nevertheless, in recent researches, utilization of smart grid

https://doi.org/10.1016/j.apenergy.2017.11.067 Received 20 August 2017; Received in revised form 25 October 2017; Accepted 11 November 2017

⁎ Corresponding author. E-mail address: [email protected] (S. Mousavizadeh).

Applied Energy 211 (2018) 443–460

Available online 22 November 2017 0306-2619/ © 2017 Elsevier Ltd. All rights reserved.

T

technologies has been propounded as an effective tool to improve the resiliency of distribution networks [12,13]. In [14] the plan develop- ment process and the R&D areas requirement in this field has been investigated. In [15] a conceptual approach for the transition process to achieve a resilient distribution system through the implementation of microgrids (MGs) has been proposed. Ref. [16] focuses on the role of networked MGs as distributed systems for enhancing the system re- siliency against extreme events. In [17] hardening strategies and smart grid technologies as tools to increase the system resiliency are dis- cussed. In [18], a decision framework is developed for grid moder- nization and the system performance improvement in restoration pro- cess under extreme weather conditions.

In addition, recent studies address the issue of quantifying the

resiliency of distribution network [19,20]. In [19] a new method have been presented for measuring the resiliency of a distribution system; the objective of this paper is to enable resiliency by comparing some pos- sible existing reconfiguration schemes and choosing the one, which increases the resiliency most. In [20], resiliency evaluation of dis- tribution system has been modeled as a multi-criteria decision making problem and calculated using graph theoretic methods.

However, the methods presented in the aforementioned references mainly discuss the topological characteristics of the distribution net- works and cannot adequately represent electrical features, existing capabilities, and the nature of the events in this part of the power grid.

Ref. [21] proposes a service restoration process that uses MGs with inadequate generation capacity to serve critical loads on distribution

Nomenclature

Indices and sets

N set of the nodes in the network Λ set of the lines in the network N set of the loads with control capabilities Ω set of the wind power uncertainties K set of the formable microgrids M set of the all distributed generations E set of the energy storages NADG set of the nodes connected to all DGs NCDG set of the nodes connected to master DGs Ψℓ set of the nodes belongs to the line ℓ

−Λi ζ x set of the lines belongs to the path x between node i and

node ζ

− ′

Λi μ z set of the lines belongs to the path ′z between node i and

node μ −Λi μ

z set of the lines belongs to the path z between node i and node μ

−AllPathi ζ set of the paths between node i and node ζ WTur set of the wind turbines NWTur set of the nodes connected to wind turbines NES set of the nodes connected to storages n index of the wind turbines i index of the buses k index of the microgrids m index of the distributed generations e index of the energy storages θL i( ) set of the connected lines to the node i

Variables and parameters

NCDGs number of the master DGs NMGs maximum number of formable MGs Pi d t

b , , blocks of the load control

PLoadi t, predicted load i at the period t ηe

ES efficiency of the storage rℓ the resistance of the line ℓ Vi k t ω, , , voltage magnitude of the nodes BigM a big number SOCe

ES initial, initial SOC of storage ζ lowest common ancestors of the parents of the node i αi k, the binary variable indicating whether node i belongs to

MG k βℓ the binary variable indicating whether the line ℓ is active HLℓ the binary variable indicating the health of the line ℓ after

the event HBi the binary variable indicating the health of the bus i after

the event

Secℓ the binary variable indicating the existence of switch on the line ℓ

−Pathi ζ x the binary variable indicating whether the path x between

node i and node ζ is active −NumPathi ζ number of the paths between node i and node ζ

P Q Q, ,m DG Max

m DG Max

m DG Min, , , the active and reactive power limits of DGs

P P Q, ,m k t DG s

m k t ω DG dep

m k t ω DG dep

, , ,

, , , ,

, , , , scheduled and generated active and reactive

power of DGs SOC SOC,e

ES Max e ES Min, , maximum and minimum SOC limits of the

storages Rate Rate,e

ch Max e dch Max, , maximum charge and discharge rate of the

storages P P Q, ,i k t

L s i k t ω L dep

i k t ω L dep

, , ,

, , , ,

, , , , scheduled and supplied active and reactive power

of the loads P P,i k t

LC s i k t ω LC dep

, , ,

, , , , scheduled and deployed active load control of the loads

σ σ,i d t L s

i d t ω L dep

, , ,

, , , , scheduled and deployed block for the load controls

PES PES,e k t dch s

e k t ω dch dep

, , ,

, , , , scheduled and deployed discharging power of the

storage PES PES,e k t

ch s e k t ω ch dep

, , ,

, , , , scheduled and deployed charging power of the

storage SOC SOC,e k t

ES s e k t ω ES dep

, , ,

, , , , scheduled and realized SOC of the storage

γ γ,e k t ES s

e k t ω ES dep

, , ,

, , , , scheduled and realized status of the storages

γn k t ω Wnd dep , , ,

, binary variable indicating whether wind turbine n is connected to MG k

PProd ,n t ω Wnd

n t Wnd Max

, , , , predicted and maximum output power of the wind

turbines P P,n k t ω

Wnd dep n k t Wnd s

, , , ,

, , , scheduled and realized active power generation of the

wind turbines PΔ i k t ω

LC , , , difference in the scheduled and deployed load control

programs PΔ m k t ω

DG , , , difference in the scheduled and generated active power of

DGs flow flow,t ω

P t ω

Q ℓ, , ℓ, , the active and reactive power flows of the line ℓ

flow flow,P Max Q Maxℓ ,

ℓ , maximum active and reactive power flow limits

of the line ℓ Vm

DG set, given voltage magnitude of the nodes connected to master DGs

Pi t LC Max ,

, maximum limit of the load control options Zflow Zflow,t ω

P t ω

Q ℓ, , ℓ, , slack variables to make the power equality con-

straints valid F1 ,F2ℓ ℓ indicating electrical characteristic of the line ℓ

−NLinei μ z number of the lines in the path −Λi μ

z

− ′

NLinei μ z number of the lines in the path −

′ Λi μ

z

−NLinei ζ x number of the lines in the −Pathi ζ

x

D number of the load control blocks Cape

ES capacity of the storage xℓ the reactance of the line ℓ δi k t ω, , , voltage angles of the nodes δ Max maximum limit of the voltage angles

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

444

feeders after an extreme event; determining the availability of the MGs for critical load restoration is the main contribution of this study. Ref. [22], investigates the optimization problem of the back-feed service restoration in smart distribution networks and optimize droop settings of DGs in islanded mode to maximize restored load.

In [23–25], optimal scheduling problem has been investigated for energy management in MGs. The objective functions of these methods is minimization of the net operational cost by considering the un- certainties of the renewable energy resources; however, the possibility of the severe event occurrences is not taken into consideration.

In [26–29], effective approaches for resilient scheduling of the static MGs are presented. However, the topological structure of the distribu- tion network and its limitations are rarely examined.

In [30,31] the subject of MG construction in active distribution networks from the perspective of energy management and optimal use of renewable energy resources is discussed. In [32], a new networked MGs structure is proposed to facilitate the MGs optimal scheduling, and the time-based demand response programs are exploited to mitigate the costs of the consumers. A novel framework for day-ahead scheduling of the smart distribution network, based on robust optimization approach, for immunizing system performance against the worst-case realization of the uncertain variables is presented in [33]. However, occurrence of the severe events is neglected in abovementioned studies.

In [34] a linear model is proposed for MGs formation in distribution networks after the disasters. It considers predetermined and fixed number of the formable MGs, and utilization of different configurations is not considered. Ref. [35], developed the algorithm presented in [34] and the impact of reconfiguration capabilities on the load restoration process is examined. The proposed model is formulated based on second-order integer programming; however, again the number of MGs and DGs as the master control unit is predetermined and fixed.

In addition, the structural and technical constraints of the dis- tribution networks, such as the limitations in the number of the sec- tionalizers for changing the status of the lines, is not properly modeled in the literature. These items can put many challenges in designing and providing the necessary measures to improve the system performance in critical conditions. Additionally, the presented formulations in [34,35] lack the failure period and the time related characteristics of the events. Considering the event-based nature of the resiliency concept and regarding the significant impacts of energy storage (ES) units on supplying the sensitive loads in critical conditions, modeling of these technologies will be beneficial in resiliency analysis.

On the other hand, due to the stochastic nature and low controll- ability of the renewable energy resources, providing a framework to scrutinize the influences of increasing the penetration level of these uncertain resources from the resiliency point of view will be of great importance. However, the literature in this realm have not addressed this issue. Moreover, owing to the significant role of the load controls and demand side programs in effective management of the system, modeling of these resources and investigating the constraints of their utilizations in critical situations is also indispensable.

Therefore, in this paper, first, a novel procedure for analyzing the resiliency of distribution networks is defined, then in order to evaluate the performance of the system in extreme events, a two-stage frame- work based on mixed-integer linear programming is presented. The proposed method attempts to provide a comprehensive formulation to resolve shortcomings of the previous researches, and maximize the system resiliency through the formation of MGs and optimal manage- ment of distributed energy resources (DERs). In the proposed approach, the topological and electrical characteristics of the distribution network are modeled as an integrated linear optimization problem. In summary, the main contributions of this paper are as follows:

• Investigating the corresponding features of the resiliency concept in distribution networks, a comprehensive formulation is proposed for analyzing and evaluating its dimensions, considering different

capabilities and constraints. In order to better diagnose the system weak points and deficiencies, the proposed framework distinguishes among the resistancy, recovery and resiliency characteristics in distribution networks, and presents separate indices for them. Employing this approach, it can be conceivable to design and utilize appropriate measures to enhance system performance in each fea- ture.

• In this paper, a linear path-based method is proposed for proper modeling of the topological characteristics of the distribution net- work, and the mentioned shortcomings in [20,34,35] are relieved. In proposed method, master control capabilities of the DGs are properly formulated and participation of these kind of DGs in MGs formation are considered within the optimization problem.

• Since resiliency is an event-based concept and severity of the events is also effective on the failure periods after the occurrence, their time related characteristics are also considered in the suggested formulation. In other words, in this paper, intensity of the events are modeled through two alternatives: the affected elements and the failure duration after the disaster. In addition, considering the time dimension, the impacts of energy-limited resources like ES units on system resiliency is also analyzed.

• Proposed formulation considers load controls, and demand side capabilities in MGs formation process and employs these opportu- nities in effective management of the system in critical situations. In addition, load curtailment through applying the demand side pro- grams is utilized to avoid complete blackout of the MGs.

• Providing a two-stage framework, stochastic nature and limited controllability of the renewable energy resources are modeled through multiple scenarios within the optimization problem, and the influence of the increasing penetration level of these uncertain resources is investigated from the resiliency concept’s perspective.

The rest of this paper is organized as follows: in Section 2, the concept of resiliency in distribution networks is explained. The pro- posed procedure for resiliency analysis and its various steps are de- scribed in Section 3. Section 4, provides the formulation of the proposed method in form of a linear optimization problem. In Section 5, several simulations have been performed to verify the validity of the proposed algorithm. Finally, concluding remarks are given in Section 6.

2. Concept of resiliency in distribution networks

Several definitions have been provided for the resilience concept over the last decade [36,37]. According to [38], “system resilience is the ability to prepare for and adapt to changing conditions, withstand, and recover rapidly from disruptions”. This definition incorporates all the active and inactive prospects of this concept (see Fig. 1).

However, this definition takes into account the resiliency

Fig. 1. Resilience features.

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

445

Fig. 2. The flowchart of the proposed two-stage method.

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

446

characteristics of the power systems in a general form. In other words, depending on the implementation context, the system under study (transmission or distribution level), as well as the horizons of the stu- dies (short-term or long-term), the way of evaluation and the approach will be different and some of the mentioned features and dimensions of this concept will be more emphasized.

For example, in transmission networks due to the interconnected structure, the system's withstanding and tolerance are more accen- tuated than its recovery feature. In other words, in transmission net- works, system ability to prevent propagation of the fault impacts is the main objective of the resiliency studies. Implementing measures like islanding of the network to avoid the cascading phenomena [39] is an example of such sort of studies.

Due to special characteristic of distribution networks, the resiliency in these systems are extremely affected by the smart grid technologies. For example, in a distribution network that includes no MGs in normal operation mode and form MGs only after the catastrophes, the recovery feature of this concept will be more prominent. In fact, in such kind of systems, due to the simple and radial structure of the grid and the heterogeneous geographical scattering of the consumers with different priorities, the ability to prevent the propagation of the failure is pos- sible only with expense of the exorbitant costs. On the other hand, in an active distribution network that employs multiple MGs in the grid- connected mode, power balancing within the MGs, at the time of con- tingency, will be realized by appropriate applying of the pre-event ac- tions such as charging the storages or reducing the loads by demand side programs. In other words, achieving the resistancy feature in dis- tribution networks would be feasible by exploiting MGs.

3. Resiliency analysis framework

Based on the previous section, a suitable procedure for resiliency analysis in distribution networks is presented here. This process consists of four steps, as follows:

• Step 1: Exploring the events As mentioned earlier, resilience is an event-based concept. In this

perspective, regarding the different nature of the events, first, based on the characteristics of the system under study, including spatial and geographic locations of the components, operation principles, opera- tor’s concerns, and existing capabilities, it is necessary to diagnose the network condition. In general, the events considered in the resilience studies include various types of natural disasters such as floods, storms, earthquakes, and human-made sabotages like a terrorist attack, vand- alism and so on.

• Step 2: Extracting the failure scenarios At this step, depending on the type and the nature of the identified

events, the impacts of these incidents will be mapped to the functioning of the network equipment. In recent years, extensive researches have been carried out on the development of the methods to analyze the associated issues. Some of these methods apply the fragility curves [40,41]; using these curves developed for each equipment, the failure probability of the elements is calculated based on the severity of the event, then comparing with the threshold value, the worst-case scenario is extracted.

• Step 3: Evaluating the system performance At this step, the aim is to investigate the behavior and response of

the system in the extracted scenarios and analyzing the different di- mensions of the resiliency. As mentioned earlier, in distribution net- works, the penetration level of smart grid technologies, and in parti- cular the MGs facility, affects the approach of resiliency assessment. In

this paper, the potential of MGs formation after the events in im- provement of the network resiliency is considered, and a novel linear model is proposed to accommodate various flexibilities.

• Step 4: Calculating the features Based on the evaluations carried out in the previous steps, the de-

sired indices will be calculated in this step to quantify the resiliency and its features in distribution network.

4. Problem formulation

The proposed method, integrates the topological and electrical characteristics of the distribution network into one linear optimization problem. Moreover, two-stage stochastic programming approach is utilized to incorporate the random nature of the renewable energy re- sources. In the first stage, the MGs formation constraints and the op- timal scheduling of the DERs are taken into account. First stage vari- ables include the topological variables (representing the status of the network buses, lines and deployment of various configurations), and electrical variables (representing the scheduled values for the DERs). The second stage characterizes the real-time operation of the MGs based on the realized uncertainties. The flowchart of the proposed method is depicted in Fig. 2. It should be noted that in this paper a scenario-based approach is employed to model the uncertainty of the output power of wind turbines. For this purpose, utilizing the historical data or applying the random generation process based on the Weibull probability density function (PDF), wind speed for the study horizon will be predicted. Then, utilizing the scenario reduction methods, the number of scenarios are reduced to diminish the computational burden. Finally, applying the power-velocity curve of the wind turbines, the output powers of the wind resources for considered scenarios will be determined. More de- tails of this modeling approach are given in [42]. A similar scenario- based model can also be applied for representing related uncertainties of the load consumptions during the failure period.

4.1. Topological constraints

In order to properly model the topological characteristics of the system, a path-based algorithm based on graph related theories is proposed. Switching operations, control capabilities of DGs and de- ployment of various configurations are also involved in this part of the model.

For this purpose, first, the network graph considering all buses and lines is extracted. Considering the tie lines, this graph will compose of several loops. The set of the nodes (buses) and the edges (lines) of the graph are represented by N and Λ, respectively. In addition, the for- mation of MGs depends on the existence of DGs with the ability to control the frequency and bus voltages. In other words, in each MG, there should be at least one DG with such capability. If there were several units with this capability in a MG, only one of these units will be chosen as the master unit. Set of the connected nodes to all DGs is represented by NADG, and NCDG is the set of the nodes connected to the units with master control capabilities.

The k-th member of NCDG is correspondent to k-th MG. In other words, in MG k, the member k of NCDG is chosen as the master unit. In addition, selecting this unit as the master, connected node to this DG, is called as the root node. For example, the sample network shown in Fig. 3 has three units with master control capability, so NCDG includes three members b m r{ , , }. Since in this network, units 1 and 3 are chosen as the master, therefore only two MGs (MG1 and MG3) are formed and the nodes b and r are considered as the root nodes. The DG unit con- nected to the node m is considered as a slave unit.

Moreover, in each MG, there should be sufficient capacity to supply the loads and achieve the power balance. If there are not enough ca- pacity or sufficient numbers of switches to form MGs with such

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

447

condition, the consumers in that sector will not be supplied. For ex- ample, as shown in Fig. 3, nodes n, a, and h do not belong to any MGs. Based on the above explanation, topological constraints can be cate- gorized as follows:

(a) Connection limit

Due to the technical constraints, each node will only be able to belong to one of the mentioned MGs or none of them. Therefore, the connection constraint of each node can be modeled as Eq. (1).

∑ ⩽ ∀ ∈ =

α i1 N k

N

i k 1

,

MGs

(1)

If the variable αi k, is equal to one, it means that node i is belong to MG k, and member k of the set NCDG is chosen as the master unit.

(b) Status of the root nodes

The hypothetical node i can be connected to MG k if only k-th member of the set NCDG (for example, node r) is chosen as the root node. This constraint can be stated as follows:

⩽ = ∀ ∈ ∀ ∈α α r k i kN ( ) , N, K.i k r k, , CDG (2)

(c) Connection of the parent nodes

To connect a node (for example i) to a MG, at least one of its parent nodes should belong to that MG. The parent nodes are the first node in the paths between the node i and the root node corresponding to that MG (for instance, the node k). For example, in Fig. 4, nodes m, n, p, and q are parents of the node i. Since none of the paths contains node a, then this node will not be construed as the parent of the node i to reach the node k. To find the parent nodes, BFS search algorithm [43] is applied. This constraint is formulated as follows.

∑⩽ ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈

α α i j ξ kN, , Ki k j ξ

j k i k, , , i k, (3)

(d) Connectivity and radiality constraints

To ensure these constraints in all MGs, only one active path should exist between each node in a MG and the corresponding root node of that MG. In order to improve the efficiency and increase the comput- ability of the proposed algorithm, idea of the first common node is introduced here. In this approach, first, the connectivity and radiality constraints are guaranteed between the node and its corresponding common node, and then these constraints are applied to the paths be- tween common node and the root node. For example, in Fig. 4, there should be only one active path to connect node i to the root node k.

Based on the common node idea, first, the connectivity and radiality constraints are checked for the paths between nodes i and ζ , and then these constraints are applied to the paths between nodes ζ and k.

To find this common node, each path between nodes i and k should be extracted. This will increase the computation burden of the proposed algorithm. To improve the efficiency, well-organized methods, devel- oped to extract the lowest common ancestor (LCA) in a graph, is uti- lized. In graph theory, the LCA of two nodes x and y in a graph is the lowest node that has both x and y as descendants [43]. For example, in Fig. 4, the node r is the LCA node of the nodes q and p. In addition, in this graph, node ζ is the LCA of the parents of the node i. To find the LCA, the method presented in [44] is also used. In this Ref., finding the LCA has been decreased to a range minimum query (RMQ) problem. Solutions to RMQ can be used to create fast answers to LCA. In order to guarantee the connectivity constraint between the hypothetical nodes i and its LCA (node ζ ), Eqs. (4) and (5) would be applied.

⩽ = ∀ ∈α α ζ ξ kLCA ( ) , Ki k ζ k i k, , , (4)

∑ ⩾ ∀ ∈ ∀ ∈ =

− −

−

α k xPath K, AllPath x

i ζ x

i k i ζ 1

NumPath

,

i ζ

(5)

The status of each path, according to (6), is related to the status of the lines in that path. The Eqs. (7) and (8) are applied for linearization of the (6).

∏= ∀ ∈− ∈

−

−

β xPath AllPathi ζ x

i ζ ℓ Λ

ℓ

i ζ x (6)

⩽ ∀ ∈ ∀ ∈− − −β xPath AllPath , ℓ Λi ζ x

i ζ i k x

ℓ (7)

∑⩾ − − ∀ ∈− ∀ ∈

− −

−

β NLine xPath ( 1) AllPathi ζ x

i ζ x

i ζ ℓ Λ

ℓ

i ζ x (8)

If there are several paths between the nodes i and ζ , the radiality constraint should be ensured as well. For this purpose, each two sets of combinations of the paths should be examined and appended the non- formation loop constraints. Similar to previous section, the idea of the first common node is applied again for each two paths. For example, in Fig. 4, the paths z and ′z , include nodes i m μ s ζ{ , , , , } and i n μ t ζ{ , , , , }, re- spectively. First, the lowest common node, between the two paths (node μ), will be identified, then the Eq. (9) is used to avoid the loop for- mation. In other words, checking the loops is continued to node μ, and the formation of the loop after this node will be checked when in- vestigating the paths between nodes μ and ζ . Eq. (9) states that in order to avoid loops, at least one of the variables indicating the statues of the corresponding lines must be set to zero.

Fig. 3. Service area of each MG.

Fig. 4. Sample graph.

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

448

∑ ∑+ ⩽ + − ∀ ∈ ∀ ∈ ∀ ′∈

′ − − ′

−

− − ′

β β NLine NLine z( 1) AllPathi μ z

i μ z

i μ ℓ Λ

ℓ ℓ Λ

ℓ i μ z

i μ z

(9)

(e) Boundary lines constraints

If the nodes of two sides of a line are not belonged to same MG, according to (10), the binary variable, indicating the status of the line, should be set to zero.

< ⊙ ∀ ∈ ∀ ∈ ∀ ∈β α α i j kℓ Λ , , Ψ , Ki k j kℓ , , ℓ (10)

In this equation, the symbol ⊙ represents the logical XNOR opera- tion between the two binary variables αi k, and αj k, . By defining the auxiliary binary variables χ kℓ, and ′χ kℓ, , the Eq. (10) could be linearize using Eqs. (11)–(13).

< + ′ ∀ ∈ ∀ ∈β χ χ kℓ Λ , Kk kℓ ℓ, ℓ, (11)

⎧

⎨ ⎩

⩽ ⩽ ⩾ + −

∀ ∈ ∀ ∈ ∀ ∈ χ α χ α χ α α

i j k 1

, Ψ , ℓ Λ , K k i k

k j k

k i k j k

ℓ, ,

ℓ, ,

ℓ, , ,

ℓ

(12)

⎧

⎨ ⎪

⎩ ⎪

′ ⩽ − ′ ⩽ − ′ ⩾ − −

∀ ∈ ∀ ∈ ∀ ∈

χ α

χ α

χ α α i j k

1

1

1 , Ψ , ℓ Λ , K

k i k

k j k

k i k j k

ℓ, ,

ℓ, ,

ℓ, , ,

ℓ

(13)

(f) Constraints of the nodes and lines status

Due to the economic issues, it is not practical to install sectionalizers or remote control switches on all lines of distribution networks. As a result, it is necessary to model the non-existence of the switch to change the status of them. On the other hand, in case of damage to the lines, inactive state should be opted for them. Therefore, Eqs. (14) and (15) model the status of a line. Eq. (16) is employed to model the failure state of the damaged nodes.

⩽ ∀ ∈β HL ℓ Λℓ ℓ (14)

⩾ − × ∀ ∈β Sec HL(1 ) ℓ Λℓ ℓ ℓ (15)

⩽ ∀ ∈ ∀ ∈α HB i kN, Ki k i, (16)

4.2. Electrical constraints

The equations for modeling the electrical constraints of the problem categorized as follows.

(a) Load constraints

In this paper, it is assumed that the system operator by utilizing small solar panels, remote control switches, or demand response con- tracts is able to directly control some of the loads. Eq. (17.1) indicates the scheduled active power at the first stage for each load considering the control capabilities. Eq. (17.2) shows the amount of active and re- active power supplied in each scenario at the second stage. Eq. (17.2) also states if the load control is applied to the active power, in ac- cordance with the power factor, the reactive power of that load is re- duced as well.

= × − ∀ ∈ ∀ ∈P α P k iPLoad K, Ni k t L s

i k i t i k t LC s

, , ,

, , , , ,

(17.1)

⎧

⎨ ⎪

⎩ ⎪

= − ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈

= × − × ∀ ∈ ∀ ∈ ∀

∈ ∀ ∈

P P P k i t ω

Q α φ P k i t

ω

Δ K, N, T, Ω

QLoad tan ( ) K, N,

T, Ω

i k t ω L dep

i k t L s

i k t ω LC

i k t ω L dep

i k i t i i k t ω LC dep

, , , ,

, , ,

, , ,

, , , ,

, , , , , ,

(17.2)

Eqs. (18.1) and (18.2) model the control capabilities of the loads at the first and second stages, respectively.

∑ ∑

∑

⎧

⎨

⎪ ⎪⎪

⎩

⎪ ⎪ ⎪

⩽ × ∀ ∈ ∀ ∈ ∀ ∈

= × ∀ ∈ ∀ ∈

⩽ ∀ ∈ ∀ ∈

∀ ∈ =

=

P α P k i t

PL σ P i t

σ i t

K, N , T

N , T

1 N , T

i k t LC s

i k i t LC Max

k i k t LC s

d

D

i d t L s

i d t b

d

D

i d t L s

, , ,

, , ,

K , ,

,

1 , ,

, , ,

1 , ,

,

(18.1)

∑

∑

⎧

⎨

⎪ ⎪⎪

⎩

⎪ ⎪ ⎪

⩽ × ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈

= × ∀ ∈ ∀ ∈ ∀ ∈

⩽ ∀ ∈ ∀ ∈ ∀ ∈

∀ ∈ =

=

P α DR k i t ω

sum P σ P i t ω

σ i t ω

K, N, T, Ω

N , T, Ω

1 N , T, Ω

i k t ω LC dep

i k i t L Max

k i k t ω LC dep

d

D

i d t ω L dep

i d t b

d

D

i d t ω L dep

, , , ,

, , ,

K , , , ,

1 , , ,

, , ,

1 , , ,

,

(18.2)

Eq. (19) shows the difference in the scheduled and deployed amount of the load control programs.

= − ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈P P P k i t ωΔ K, N , T, Ωi k t ω LC

i k t ω LC dep

i k t LC s

, , , , , , ,

, , ,

(19)

(b) DGs constraints

Eqs. (20.1) and (20.2) represent the active and reactive power limits of DGs in the first and second stages, respectively.

⩽ × ∀ ∈ ∀ ∈ ∀ ∈ =P α P m k t r mM, K, T, N ( )m k t DG s

r k m DG Max

, , ,

, ,

ADG

(20.1)

⎧

⎨

⎪ ⎪ ⎪ ⎪⎪

⎩

⎪ ⎪ ⎪ ⎪ ⎪

⩽ × = ∀ ∈ ∀ ∈ ∀

∈ ∀ ∈ = + ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈

⩽ × = ∀ ∈ ∀ ∈ ∀

∈ ∀ ∈ ⩾ × = ∀ ∈ ∀ ∈ ∀

∈ ∀ ∈

P α P r m m k t

ω P P P m k t ω

Q α Q r m m k t

ω Q α Q r m m k t

ω

N ( ) , M, K,

T, Ω Δ M, K, T, Ω

N ( ) , M, K,

T, Ω N ( ) , M, K,

T, Ω

m k t ω DG dep

r k m DG Max

DG

m k t ω DG dep

m k t DG s

m k t ω DG

m k t ω DG dep

r k m DG Max

DG

m k t ω DG dep

r k m DG Min

DG

, , , ,

, ,

, , , ,

, , ,

, , ,

, , , ,

, ,

, , , ,

, ,

(20.2)

(c) Energy storages constraints

The Eqs. (21.1), (21.2) show the charging and discharging limits of the storages in the first and second stages, respectively.

⎧

⎨

⎪ ⎪

⎩ ⎪ ⎪

⩽ = ∀ ∈ ∀ ∈ ∀ ∈

⩽ × ∀ ∈ ∀ ∈ ∀ ∈

⩽ − × ∀ ∈ ∀ ∈ ∀ ∈

⩽ × ∀ ∈ ∀ ∈ ∀ ∈ =

γ α i e e k t

PES γ Rate e k t

PES γ Rate e k t

PES α Rate e k t i e

N ( ) , E, K, T

E, K, T

(1 ) E, K, T

E, K, T, N ( )

e k t ES s

i k

e k t ch s

e k t ES s

e ch Max

e k t dch s

e k t ES s

e dch Max

e k t dch s

i k e dch Max

, , ,

, ES

, , ,

, , , ,

, , ,

, , , ,

, , ,

, ,

ES

(21.1)

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

449

⎧

⎨

⎪ ⎪ ⎪ ⎪

⎩

⎪ ⎪ ⎪ ⎪

⩽ = ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈

⩽ × ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈

⩽ − × ∀ ∈ ∀ ∈ ∀

∈ ∀ ∈ ⩽ × = ∀ ∈ ∀ ∈ ∀

∈ ∀ ∈

γ α i e e k t ω

PES γ Rate e k t ω

PES γ Rate e k t

ω PES α Rate i e e k t

ω

N ( ) , E, K, T, Ω

E, K, T, Ω

(1 ) E, K,

T, Ω N ( ) , E, K,

T, Ω

e k t ω ES dep

i k

e k t ω ch dep

e k t ω ES dep

e ch Max

e k t ω dch dep

e k t ω ES dep

e dch Max

e k t ω dch dep

i k e dch Max

, , , ,

, ES

, , , ,

, , , , ,

, , , ,

, , , , ,

, , , ,

, ,

ES

(21.2)

Eqs. (22.1) and (22.2) also indicate restrictions of the maximum, minimum and initial charge levels of the storages, respectively.

⎧

⎨ ⎪

⎩ ⎪

⩽ × = ∀ ∈ ∀ ∈ ∀ ∈

⩾ × = ∀ ∈ ∀ ∈ ∀ ∈

= ∀ ∈ ∀ ∈=

SOC α SOC i e e k t

SOC α SOC i e e k t

SOC SOC e k

N ( ) , E, K, T

N ( ) , E, K, T

| E, K

e k t ES s

i k e ES Max

e k t ES s

i k e ES Min

e k t ES s

t e ES initial

, , ,

, ,

ES

, , ,

, ,

ES

, , ,

0 ,

(22.1)

⎧

⎨

⎪ ⎪ ⎪

⎩

⎪ ⎪ ⎪

⩽ × = ∀ ∈ ∀ ∈ ∀

∈ ∀ ∈ ⩾ × = ∀ ∈ ∀ ∈ ∀

∈ ∀ ∈ = ∀ ∈ ∀ ∈ ∀ ∈=

SOC α SOC i e e k t

ω SOC α SOC i e e k t

ω SOC SOC e k ω

N ( ) , E, K,

T, Ω N ( ) , E, K,

T, Ω | E, K, Ω

e k t ω ES dep

i k e k t ES Max

e k t ω ES dep

i k e k t ES Min

e k t ω ES dep

t e ES initial

, , , ,

, , , ,

ES

, , , ,

, , , ,

ES

, , , ,

0 ,

(22.2)

Eqs. (23) and (24) compute the SOC of the storages in various time intervals in the first and second stages.

= + × − ×

= ∀ ∈ ∀ ∈ ∀ ∈

+ × SOC SOC PES PES

i e e k tN ( ) , E, K, T

e k t ES s

e k t ES s

e k t ch s η

Cap e k t dch s

η Cap, , 1 ,

, , ,

, , ,

, , , 1

ES

e ES

e ES

e ES

e ES

(23)

= + × − ×

= ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈

+ × SOC SOC PES PES

i e e k t ωN ( ) , E, K, T, Ω

e k t ω ES dep

e k t ω ES dep

e k t ω ch dep η

Cap e k t ω dch dep

η Cap, , 1, ,

, , , ,

, , , ,

, , , , 1

ES

e ES

e ES

e ES

e ES

(24)

(d) Wind turbines constraints

Eqs. (25.1) and (25.2) model the constraints of the wind resources output power with respect to the predicted scenarios.

⩽ × = ∀ ∈ ∀ ∈ ∀

∈

P α P i n n k tN ( ) , WTur, K,

T n k t Wnd s

i k n t Wnd Max

WTur, , ,

, , ,

(25.1)

⎧

⎨

⎪ ⎪

⎩ ⎪ ⎪

⩽ = ∀ ∈ ∀ ∈ ∀ ∈ ∀

∈ = × ∀ ∈ ∀ ∈ ∀ ∈ ∀

∈

γ α i n n k t ω

P γ n k t ω

N ( ) , WTur, K, T,

Ω Prod WTur, K, T,

Ω

n k t ω Wnd dep

i k WTur

n k t ω Wnd dep

n k t ω Wnd dep

n t ω Wnd

, , , ,

,

, , , ,

, , , ,

, ,

(25.2)

(e) Power balance in each MG

Eq. (26) ensures the active power balance at each time interval and at all MGs in the first stage.

∑ ∑ ∑ ∑ ∑+ + − −

= ∀ ∈ ∀ ∈ = = = = =

P P PES PES P

k t0 K, T m

N

m k t DG s

n

N

n k t Wnd s

e

N

e k t dch s

e

N

e k t ch s

i

N

i k t L s

1 , ,

,

1 , ,

,

1 , ,

,

1 , ,

,

1 , ,

, DG Wtr ES ES Load

(26)

(f) Line flow limits

Eqs. (27) and (28) model the constraints of the active and reactive power flows of the lines.

− × ⩽ ⩽ × ∀ ∈ ∀ ∈ ∀

∈

β flow flow β flow t ωℓ Λ , T,

Ω

P Max t ω

P P Max ℓ ℓ

, ℓ, , ℓ ℓ

,

(27)

− × ⩽ ⩽ × ∀ ∈ ∀ ∈ ∀

∈

β flow flow β flow t ωℓ Λ , T,

Ω

Q Max t ω

Q Q Max ℓ ℓ

, ℓ, , ℓ ℓ

,

(28)

(g) Power balance in each node

Eqs. (29) and (30), respectively, state constraints of the active and reactive power balance in each node at the second stage.

∑

∑

+ + − −

=

∀ ∈ ∀ ∈ ∀ ∈

∀ ∈

∈

P P PES PES P

flow

i t ω

[ ( ) ]

N, T, Ω

k m k t ω DG dep

n k t ω Wnd dep

e k t ω dch dep

e k t ω ch dep

i k t ω L dep

θL i t ω

P

K , , ,

, , , ,

, , , ,

, , , ,

, , , ,

,

ℓ ( ) ℓ, ,

(29)

∑ ∑− = ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈ ∈

Q Q flow i t ω[ ] N, T, Ω k

m k t ω DG dep

i k t ω L dep

θL i t ω

Q

K , , ,

, , , ,

,

ℓ ( ) ℓ, ,

(30)

(h) Node voltage constraints

Eqs. (31) and (32) model the limits of the node voltage magnitude and angles, respectively.

× ⩽ ⩽ × ∀ ∈ ∀ ∈ ∀ ∈ ∀ ∈α V α i k t ω0.9 1.1 N, K, T, Ωi k i k t ω i k, , , , , (31)

− × ⩽ ⩽ × ∀ ∈ ∀ ∈ ∀ ∈ ∀

∈

α δ δ α δ i k t ωN, K, T,

Ω i k

Max i k t ω i k

Max , , , , ,

(32)

Eq. (33) also states that if the DG k is selected as a master unit, the voltage magnitude of the corresponding node =i k( N ( ) )CDG is also set at the controlled value. Eq. (34) sets the voltage angle of the mentioned node to zero.

= × ∀ ∈ ∀ ∈ ∀ ∈V α V i t ωN , T, Ωi k t ω i k k DG set

CDG, , , , , (33)

− − × ⩽ ⩽ − × ∀ ∈ ∀ ∈ ∀

∈

α δ δ α δ i t ω(1 ) (1 ) N , T,

Ω i k

Max i k t ω i k

Max CDG, , , , ,

(34)

(i) Load flow constraints

In this paper, the introduced method in [45] is used to conduct load flow calculations in distribution networks. In this Ref., an acceptable linear approximation is applied to calculate the magnitude and the angle of the bus voltages. The mentioned method also takes into ac- count the reactive power flows of the lines. According to [45], linear- ized power flow equations in distribution networks could be modeled by (35)–(37).

= +

= +

r r x

F x

r x F1 , 2ℓ

ℓ

ℓ 2

ℓ 2 ℓ

ℓ

ℓ 2

ℓ 2 (35)

∑ ∑= + − × + −

× ∀ ∈ ∀ ∈ ∀ ∈

∀ ∈ ∀ ∈ flow Zflow δ δ V V

t ω

[ ] F2 [ ]

F1 ℓ Λ , T, Ω

t ω P

t ω P

k j k t ω i k t ω

k j k t ω i k t ωℓ, , ℓ, ,

K , , , , , , ℓ

K , , , , , ,

ℓ

(36)

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

450

∑ ∑= + − × + −

× ∀ ∈ ∀ ∈ ∀ ∈

∀ ∈ ∀ ∈ flow Zflow δ δ V V

t ω

[ ] F1 [ ]

F2 ℓ Λ , T, Ω

t ω Q

t ω Q

k i k t ω j k t ω

k j k t ω i k t ωℓ, , ℓ, ,

K , , , , , , ℓ

K , , , , , ,

ℓ

(37)

Eqs. (38) and (39) also show the limits of the slack variables to make equality constraints valid when the nodes of the both sides of a line do not belong to same MG.

− − × ⩽ ⩽ − × ∀ ∈ ∀ ∈ ∀

∈

β Zflow β t ω(1 ) BigM (1 ) BigM ℓ Λ , T,

Ω t ω

P ℓ ℓ, , ℓ

(38)

− − × ⩽ ⩽ − × ∀ ∈ ∀ ∈ ∀

∈

β Zflow β t ω(1 ) BigM (1 ) BigM ℓ Λ , T,

Ω t ω

Q ℓ ℓ, , ℓ

(39)

4.3. Objective function

Since the main objective of this paper is to examine the performance of distribution networks in the face of extreme events, according to (40), the objective function of the optimization problem is maximizing total power of the restored loads with regard to their priorities.

∑ ∑ ∑ ∑= ⎛ ⎝ ⎜ ×

⎡

⎣ ⎢ − ×

⎤

⎦ ⎥

⎞

⎠ ⎟

∈ ∈ ∈ ∈

Max ObjF P π P: Pr Δ t T i

i t L

k i k t L s

ω ω i k t ω

LC

N ,

K , ,

,

Ω , , ,

(40)

The system resistancy (as one of the features of the resiliency con- cept), is the ratio of the summation of active powers of the non-inter- rupted costumers to the total active powers of the system loads, taking into account the priority factor. According to Eq. (41), this feature states the ability of the system to tolerate the event and prevent pro- pagation of its effects.

= ∑ ×

∑ × ∈

∈

−Resistancy (Pr PLoad )

(Pr PLoad ) i i t

L i t

i i t L

i t

N , ,

N , ,

Not Intrrupted

(41)

As mentioned in the Section 3, if there were connected MGs in distribution network, by implementing the appropriate preparation, it could be possible to reduce the number of the affected loads and as a result, increase resistancy of the system.

In addition, recovery index could be extracted using Eq. (42), and calculating the ratio of the expected recovered energy of the interrupted loads, to sum of the total energy of the interrupted loads in the study period.

Fig. 5. Single line diagram of the test network 1.

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

451

= ∑ ∑ ×∈ ∈

ObjF Recovery

(Benf PLoad )t T i i t L

i tN , ,Interrupted (42)

Resiliency of the system is also calculated using Eq. (43), and is the ratio of the expected supplied energy of the loads during study period (loads connected to non-faulted feeders as well as loads recovered through forming MGs) to total energy of the system loads (whether interrupted or non-interrupted).

R = + ∑ ∑ ×

∑ ∑ × ∈ ∈

∈ ∈

−esilience ObjF (Benf PLoad )

(Benf PLoad ) t T i i t

L i t

t T i i t L

i t

N , ,

N , ,

Not Intrrupted

(43)

As mentioned earlier, due to the severity of the events considered in resiliency studies, and weak structure of distribution networks, usually NInterrupted includes large portion of the loads and hence, the concept of the resiliency in distribution networks, is mainly indicative of the

recovery capabilities.

5. Numerical study

In this section, performance of the proposed algorithm for resiliency analysis in distribution networks is investigated. The simulations have been executed on a PC with Intel Core i7 CPU @4 GHz and 16 GB RAM. The proposed model has been solved using IBM ILOG CPLEX 12.6 [46] and MATLAB software, as well as considering a mip gap of 0.01%.

5.1. Test network 1

A modified 11-kV distribution network as shown in Fig. 5, com- posed of 3 feeders, 118 buses, 3 breakers, 30 sectionalizers, and 9 tie lines, is considered as the test network 1. The total active and reactive power consumptions of this network are 22.71 MW and 17.04 MVAr, respectively, and the uncertainty of the loads is neglected. The detailed data about the system peak loads and line parameters are given in [47].

Required information about the location and the capacity of DGs, wind turbines and energy storages is given in Table 1. In this network, 7 units (DG1 to DG7) of the 14 existing DGs, have the ability to perform as a master unit, and therefore, the maximum number of the formable MGs, is equal to 7. The capacity, charging and discharging rates, effi- ciency, and initial SOC of the ES units are assumed to be 500 kWh, 100 kW, 85% and 60%, respectively. The capacity of the wind turbines is also considered 500 kW. The information of wind power uncertainties has been extracted from [29] and is as shown in Fig. 6. The daily load profile based on the P.U. value is shown in Fig. 7, and priority of the supplying loads is set to one. The duration of the switching time is also considered one hour.

Table 1 DGs, wind turbines and ES required information.

Unit Bus Cap. (kW) Qmax (kVAr) Qmin (kVAr) Unit Bus Unit Bus

DG1 17 800 600 −600 ES1 12 WT1 14 DG2 24 1000 800 −800 ES2 26 WT2 27 DG3 51 900 800 −800 ES3 40 WT3 42 DG4 59 1200 1000 −1000 ES4 53 WT4 53 DG5 67 1000 800 −800 ES5 72 WT5 74 DG6 76 1100 800 −800 ES6 80 WT6 84 DG7 107 1500 1200 −1200 ES7 90 WT7 88 DG8 7 200 150 −150 ES8 110 WT8 96 DG9 33 300 200 −200 WT9 99 DG10 43 300 200 −200 WT10 112 DG11 88 200 100 −100 DG12 103 500 300 −300 DG13 113 300 150 −150 DG14 117 300 150 −150

Fig. 6. Output power scenarios of wind resources.

Fig. 7. Daily load profile.

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

452

5.1.1. Studied events In this paper, occurrence of the flood and a vandalism have been

investigated. In first event (flood), it is assumed that by applying fra- gility curves for the power distribution poles, buses 2, 3, 4, 10, 63, 64, 65 and 89 will be in failure state. In this event, connections of the feeders 1 and 2 are interrupted from the main substation and only feeder 3 will be connected to the upstream network. In the second event (destruction of the main substation), only bus number 1 will be in failure state and connection of all three feeders to the upstream network will be broken. Duration of the failure in the first event is considered to be 15 h (between 10:00 and 24:00) and in second incident it is 20 h (between 3:00 and 22:00).

5.1.2. Impact of the switching capabilities In this case, the impacts of the switching capabilities and utilization

of tie lines on system resiliency during the failure period have been examined. Table 2 shows resistancy, recovery and resiliency indices of the test network in both events.

Since in the first event, only the loads located in feeders 1 and 2 are affected, the resistancy of the system is equal to 22.23%. In other words, about 77% of the system loads have experienced interruption at the moment of the first event occurrence. In second event, given the fact that feeder 3 is also affected, the resistancy index has fallen to zero. In addition, according to Eqs. (42) and (43), and because of the total loads interruption in the second event, recovery and resiliency indices are equal. As mentioned in previous sections, recovery is the system ability to restore the loads in the affected area (feeders 1 and 2 in the first event, and the entire system in the second event). Based on the obtained results, in the first event, about 19% of the benefits of the loads can be recovered by forming the MGs in feeders 1 and 2. By utilizing the ex- isting tie lines, this index has also increased to 32%, indicating the importance of this flexibility in improving the recovery process. These values, in the second event, are 20% and 31% respectively. In Fig. 8, the formation of the MGs in feeders 1 and 2, considering the tie lines for the first event is shown.

The resiliency index states that the system under study will lose about 63% and 80% of the welfares of its loads, due to the mentioned events. Based on the results, it is obvious that by using appropriate switching operation and applying tie lines, resiliency index of the system can be improved about 10% in these events.

Additionally, the obtained execution times in Table 2, verify the computational efficiency of the proposed linear solution for resiliency analysis in large-scale distribution networks.

Table 2 Performance of the system in both events.

Event 1 Event 2

No ties With ties No ties With ties

Execution time (sec) 3.3 24.41 7.76 352.24 Number of formed MGs 4 3 5 4 Master DGs 2, 3, 4, 6 1, 3, 6 2, 3, 4, 6, 7 1, 3, 6, 7 Resistancy 0.2223 0.2223 0 0 Recovery 0.1939 0.3200 0.2001 0.3184 Resiliency 0.3731 0.4712 0.2001 0.3184

Fig. 8. MGs formation in feeders 1 and 2 for event 1.

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

453

Moreover, to evaluate the efficiency of the proposed algorithm in comparison with the previous studies, the presented methods in [34,35] are also implemented on the test network and the obtained results are shown in Figs. 9and 10. As it can be seen from the results, the presented method in this paper has a higher performance than the proposed methods in [34,35]. The reason is that in those references, the number of MGs is predetermined and each DGs is considered as a master unit; while by varying in the configuration of DGs, the determination of the optimal configuration of MGs will be affected. Therefore, presented method in previous articles do not utilize the maximum potential of the existing capabilities; but this issue has been considered in this paper and amended. As it can be seen from Fig. 8, the ability of the presented method of this paper in considering the participation of the master sources in formation of MGs, have led to form three MGs in feeders 1 and 2. In addition, loads 94, 95, and 96 are connected to feeder 3 by closing the switch on the tie line T8.

5.1.3. Impact of load control capabilities In this case, for the loads with the consumption peak higher than

400 kW (buses 20, 32, 39, 43, 44, 45, 58, 71, 74, 85, 86, 91, 101, 102, 103, 107, 109 and 113), controlling options at five levels are

considered. In first four levels, the load is reduced 100 kW in each step, and finally in the last step, the load is interrupted completely. Table 3 shows the resistancy, recovery and resiliency indices of the studied system for this case. In this table, the switching through the tie lines is also considered.

It can be inferred from the results that in this case, with less number of MGs, a large portion of loads are recovered. Regarding that the load control capabilities do not have an influence on the resistancy of the

Fig. 9. Comparison of the obtained results for the event 1.

Fig. 10. Comparison of the obtained results for the event 2.

Table 3 Performance of the system in both events considering load control capabilities.

Event 1 Event 2

Execution time (sec) 259.1 645.3 Number of formed MGs 1 1 Master DGs 1 5 Resistancy 0.2223 0 Recovery 0.6464 0.6415 Resiliency 0.7250 0.6415

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

454

system against the propagation of the failure, this index has not changed. In this case, using the load control options and a more ef- fective management, the system resiliency index in the first and the second event has increased 25 and 32 percent, respectively; this in- dicates the significant role of the demand side resources in improving system resiliency. The reason for this eminent enhancement is that in this case, using the load control capabilities, maximum utilization of the uncertain wind turbines and the other resources have been achieved.

However, by increasing the number of variables and problem con- straints in this case, the execution times have increased for both events. It should be noted that, due to the scale of the studied distribution network and simultaneous consideration of the various capabilities in the evaluation process, these execution times are quite logical and the proposed algorithm is also applicable for implementing a recovery plan for real distribution networks.

In Figs. 11 and 12, utilization manner of the load control cap- abilities during the failure period of the studied events is represented. The system operator can employ these results for appropriate devel- opment of the control facilities in the network. For example, im- plementing the control options for the loads located in buses 39, 58, 71, 85, and 103 would not be effective adequately on system performance for the studied critical conditions. On the contrary, because of the higher amount of deployed control options in the loads 74, 86, 101, 102, 107, 109, and 103, developing the controlling facility for these loads can be more beneficial through the resiliency point of view.

In addition, the average SOC levels of the ESs that are connected to MGs, for the events 1 and 2, are represented in Figs. 13 and 14, re- spectively. ES3 in both events and ES7 in the first event are not

belonged to any MGs. Based on the results, ESs are employed for sup- plying the loads and managing wind power uncertainties, and generally the average SOC levels of them are decreased to the minimum level by the end of the failure period in both events. Furthermore, at some times, especially during the peak hours of the wind power scenarios, average SOC level of the ESs has increased, and even exceeded the initial SOC level. This is more evident in Fig. 13.

In other words, in studied events, ESs are employed to absorb the generated power of the wind resources in these periods and discharge the stored energy in the next hours. Therefore, in order to maximize utilization of the renewable energy resources and to effective man- agement of their fluctuations, it is necessary to consider absorptive flexibilities of ESs. In fact, maximum initial SOC does not necessarily mean better conditions for system resiliency. This issue emphasizes the optimal determination of the initial SOC of ESs for better management of the uncertainties of renewable resources during the failure period.

5.1.4. Impact of the renewable energy resources In Table 4, values of the objective function and recovery index of the

system in the first event for various capacities of wind turbines are presented. In Fig. 15, system resiliency index is also shown. In this case, load control capabilities are considered too. Base on the obtained re- sults, by increasing the penetration level of the wind power resources, more loads are supplied during the failure period, and hence, recovery and resiliency indices have been improved. However, simulation results demonstrate that, increasing the capacity of wind turbines more than 750 kW, has not improved the indices significantly. This result em- phasizes appropriate developing rate of the control capabilities in

Fig. 11. Average load control utilization during the scheduling period of the first event.

Fig. 12. Average load control utilization during the scheduling period of the second event.

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

455

proportion to penetration level of the uncertain sources. In addition, restored loads, and recovery and resiliency indices of

the system have been decreased by adding the capacity of the wind turbines from 1250 kW to 1500 kW. The increment in the volatility and variations of the output power of these resources has led to this result. In other words, by increasing the capacity of wind turbines to 1500 kW, there would not be enough control capabilities to moderate these fluctuations. It can be concluded that, by increasing the penetration

level of these uncertain resources, more flexibilities like ESs will be required to relieve their variation and enhance resiliency.

Another result is that, if because of the weather conditions, utili- zation of the renewable energy resources would be not possible, ac- cording to Table 4, for the penetration level of 44%, the system re- covery and resiliency will decrease about 25% and 20%, respectively. In other words, increasing the dependency of the distribution networks to these unreliable sources will put the system restoration capabilities

Fig. 13. The average SOC level of the energy storages during the first event.

Fig. 14. The average SOC level of the energy storages during the second event.

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

456

and its resiliency levels at serious challenges. This emphasizes that development of renewable energy resources should also be investigated from the resiliency perspective.

5.2. Real distribution network

In this section, the proposed algorithm is implemented on a real 20- kV distribution network of Sa'adat-Abad district of Tehran, Iran. Single line diagram of this network based on GIS data is shown in Fig. 16. This network includes one 63/20 kV substation, 14 MV feeders, 492 buses, 516 lines, 79 sectionalizers and 25 tie lines. The total active and re- active power consumptions of this network are 46.739 MW and 15.063 MVAr, respectively. The detailed data about the system loads, daily consumption pattern, and the line parameters are given in [48].

In addition, 17 DG units with a total capacity of 24.1 MW, 23 wind turbines (with capacity of 500 kW), and 17 ES units are added to the network. The required information about the location and character- istics of each resource is also given in [48]. The stochastic nature of the wind power resources is modeled based on the probability distribution functions which have been obtained from the historical data of the wind speed for Tehran [49]. Then, based on the presented approach in Sec- tion 4, the output power scenarios of these resources are generated for applying in simulations.

Due to the spatial features and geographical characteristics of the district, occurrence of the storm has been investigated for this network. It is assumed that by applying fragility curves, 36 power distribution lines will be in failure state because of the storm. The duration of the failure is considered to be 13 h (between 15:00 and 03:00), and the duration of the switching time is assumed to be one hour.

Fig. 15. The obtained resilience index for different wind resource penetration.

Table 4 System performance in different wind resource penetrations.

Turbine capacity (kW)

Penetration level (%)

Objective function (kWh)

Recovery index

0 0 103,051 0.4426 250 11 127,354 0.5470 500 22 150,482 0.6464 750 33 158,423 0.6805 1000 44 162,367 0.6974 1250 55 164,424 0.7062 1500 66 159,382 0.6846

Fig. 16. GIS map of the studied real network.

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

457

5.2.1. Simulation results In this section, as shown in Table 5, impact of six different strategies

on system performance has been studied. In Fig. 17, calculated system resiliency, recovery, and resistancy indices are represented for the mentioned strategies. As it can be seen from the results, appropriate employment of the generation resources and tie lines can significantly enhance the performance of the system against the severe events. For example, simultaneous utilization of DG units and tie lines have im- proved system recovery and resiliency about 53% and 70%, respec- tively. This necessitates more attention to the placement and allocation problem of the line switching and DGs from the resiliency point of view. Similar to the previous test network, the mentioned strategies will not have an effect on the resistancy index, and therefore, it has remained unchanged for all strategies.

In addition, these results show that in order to increase the system resiliency and improve the restoration capabilities, an appropriate strategy should be chosen based on the cost-benefit analysis. For ex- ample, in this test system, by increasing the capacity of DG units to 50%, the recovery and resiliency of the system improved about 10% and 12%, respectively, while these indices increased 38% and 29% just by adding the tie lines. Furthermore, by applying the tie lines, due to the better participation of the resources in supplying the system loads, number of formed MGs, has been decreased too.

However, for more analysis, recovery and resiliency indices for each feeders of the test network are calculated. These results are represented in Tables 6 and 7. The results show that due to the limited number of the sectionalizer switches on the real distribution networks, increasing the capacity of the DG units will not necessarily improve the restoration capabilities. For example, in feeders F1, F2, F4, F8, F9, F13 and F14, 50 percent increase in DGs capacities, without considering the tie lines, did not improve the recovery indices; while by considering the tie lines due to better resources management, recoverability of the feeders has in- creased significantly.

Furthermore, in some feeders, the strategy of adding the tie lines has improved the resiliency of the system more than the strategy of the increasing DGs capacities. For example, in feeder F1, by adding the tie lines, about 92 percent of the disconnected loads are restored, while by increasing the capacity of DGs to 100 percent, due to the limited number of sectionalizers, only 50 percent of the loads are supplied during the failure period.

However, this process is reverse for feeder F10. In other words, in feeder F10, increasing the capacity of DGs has more effect in feeder’s

Table 5 System performance in six different strategies.

With tie & 100% DG

With tie & 50% DG

With tie & 0% DG

No tie & 100% DG

No tie & 50% DG

No tie & 0% DG

Objective function (kWh)

3.63E +5

2.83E +5

2.49E +5

2.84E +5

1.55E +5

1.10E +5

Solving time (min)

21.3 19.4 16.22 5.6 3.32 0.01

Resiliency 0.7763 0.6095 0.5381 0.6104 0.3413 0.2489 Resistivity 0.2523 0.2523 0.2523 0.2523 0.2523 0.2523 Recovery 0.7008 0.4777 0.3822 0.4835 0.1236 0 Number of

MGs 10 5 0 15 8 0

Fig. 17. System resiliency, recovery, resistancy indices in six different strategies.

Table 6 Recovery indices of the network feeders.

With tie & 100% DG

With tie & 50% DG

With tie & 0% DG

No tie & 100% DG

No tie & 50% DG

No tie & 0% DG

F1 0.9276 0.9276 0.9276 0.5069 0 0 F2 0.9276 0.6186 0 0.9276 0 0 F3 0.9276 0.9276 0.9276 0.9276 0.6195 0 F4 0.5869 0.2388 0.2388 0.3481 0 0 F5 0.7071 0.1985 0 0.7071 0.1985 0 F6 0.8125 0.8125 0.8125 0.6894 0.2375 0 F7 0.5239 0.2952 0 0.5239 0.2952 0 F8 0.9276 0.1018 0.1018 0.8258 0 0 F9 0.5939 0.3218 0.3218 0.272 0 0 F10 0.9276 0.5355 0 0.9276 0.5355 0 F11 0.5013 0.0828 0.0238 0.3016 0.059 0 F12 0.3029 0.2062 0.1073 0.1956 0.099 0 F13 0.9276 0.9136 0.9136 0 0 0 F14 0.7982 0.7982 0.7982 0.3645 0 0

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

458

performance than the utilization of the tie lines. Based on these results, according to the conditions of each feeder, spatial and geographical location, number of existing switches, severity of the event and number of the damaged elements, a suitable strategy for each feeder should be developed in order to improve its performance and resiliency.

6. Conclusion

This paper presented a novel two-stage stochastic framework for analyzing the resiliency of distribution networks in the face of disasters and the proposed model has been verified on two case studies. Simulation results show that the model can utilize full potentials of the smart grids facilities such as MGs and DERs to improve the resiliency and recovery of the grid. The proposed method has the ability to con- sider the participation of master DGs in formation of dynamic MGs and therefore, this issue has led to better load restoration and system re- covery in comparison with the previous methods. Results of considering the impacts of demand side programs show that implementation of load controls in some nodes can significantly help the operator to enhance system performance in recovery stage. In addition, the impacts of the renewable energies developments on the resiliency features are also investigated. The results emphasize that these uncertain resources can put system in many challenges; and thus, further examinations are re- quired for optimal growth penetration of them in electric power sys- tems.

Due to the significant growth of natural disasters in different countries in recent years, the method presented in this paper can have important consequences in operation and planning analysis of real distribution networks. Since the method distinguishes between the re- siliency features (recovery and resistancy) and calculates separate in- dices for them, identifying the system weak points from the resiliency point of view will be more plausible and beneficial.

Based on this study, the future research efforts can include the de- termination of the optimal buses for load controls, optimal placement of DGs from the resiliency point of view and estimating the required emergency budgets for responding and restoring the system from the disruptions.

References

[1] Aki H. Demand-side resiliency and electricity continuity: experiences and lessons learned in Japan. Proc IEEE 2017;105:1443–55. http://dx.doi.org/10.1109/JPROC. 2016.2633780.

[2] Gholami A, Aminifar F, Shahidehpour M. Front lines against the darkness: enhan- cing the resilience of the electricity grid through microgrid facilities. IEEE Electrification Mag 2016;4:18–24. http://dx.doi.org/10.1109/MELE.2015. 2509879.

[3] Panteli M, Mancarella P. Influence of extreme weather and climate change on the resilience of power systems: impacts and possible mitigation strategies. Electr Power

Syst Res 2015;127:259–70. [4] Panteli M, Mancarella P. The grid: stronger, bigger, smarter?: Presenting a con-

ceptual framework of power system resilience. IEEE Power Energy Mag 2015;13:58–66. http://dx.doi.org/10.1109/MPE.2015.2397334.

[5] Wang Y, Chen C, Wang J, Baldick R. Research on resilience of power systems under natural disasters – a review. IEEE Trans Power Syst 2016;31:1604–13. http://dx. doi.org/10.1109/TPWRS.2015.2429656.

[6] Panteli M, Trakas DN, Mancarella P, Hatziargyriou ND. Power systems resilience assessment: hardening and smart operational enhancement strategies. Proc IEEE 2017;105:1202–13. http://dx.doi.org/10.1109/JPROC.2017.2691357.

[7] Espinoza S, Panteli M, Mancarella P, Rudnick H. Multi-phase assessment and adaptation of power systems resilience to natural hazards. Electr Power Syst Res 2016;136:352–61.

[8] Arghandeh R. The local team: leveraging distributed resources to improve resi- lience. IEEE Power Energy Mag 2014;12:76–83. http://dx.doi.org/10.1109/MPE. 2014.2331902.

[9] Francis R, Bekera B. A metric and frameworks for resilience analysis of engineered and infrastructure systems. Reliab Eng Syst Saf 2014;121:90–103.

[10] Panteli M, Mancarella P, Trakas D, Kyriakides E, Hatziargyriou N. Metrics and quantification of operational and infrastructure resilience in power systems. IEEE Trans Power Syst 2017;99:1. http://dx.doi.org/10.1109/TPWRS.2017.2664141.

[11] Ouyang M, Osorio L. Multi-dimensional hurricane resilience assessment of electric power systems. Struct Saf 2014;48:15–24.

[12] Liu X, Shahidehpour M, Li Z, Liu X, Cao Y, Bie Z. Microgrids for enhancing the power grid resilience in extreme conditions. IEEE Trans Smart Grid 2017;8:589–97. http://dx.doi.org/10.1109/TSG.2016.2579999.

[13] Ren L, Qin Y, Li Y, Zhang P, Wang B, Luh B, Han S, Orekan T, Gong T. Enabling resilient distributed power sharing in networked microgrids through software de- fined networking. Appl Energy 2017. http://dx.doi.org/10.1016/j.apenergy.2017. 06.006. ISSN 0306-2619.

[14] Ton DT, Wang WTP. A more resilient grid: the U.S. Department of Energy joins with stakeholders in an R&D plan. IEEE Power Energy Mag 2015;13:26–34. http://dx. doi.org/10.1109/MPE.2015.2397337.

[15] Estévez G, Espinosa A, Behnke R, Lanuzza L, Velázquez N. Achieving resilience at distribution level: learning from isolated community microgrids. IEEE Power Energy Mag 2017;15:64–73. http://dx.doi.org/10.1109/MPE.2017.2662328.

[16] Li Z, Shahidehpour M, Aminifar F, Alabdulwahab A, Turki Y. Networked microgrids for enhancing the power system resilience. Proc IEEE 2017;105:1289–310. http:// dx.doi.org/10.1109/JPROC.2017.2685558.

[17] Bie Z, Lin Y, Li G, Li F. Battling the extreme: a study on the power system resilience. Proc IEEE 2017;105:1253–66.

[18] Chen C, Wang J, Ton D. Modernizing distribution system restoration to achieve grid resiliency against extreme weather events: an integrated solution. Proc IEEE 2017;105:1267–88. http://dx.doi.org/10.1109/JPROC.2017.2684780.

[19] Chanda S, Srivastava AK. Defining and enabling resiliency of electric distribution systems with multiple microgrids. IEEE Trans Smart Grid 2016;7:2859–68. http:// dx.doi.org/10.1109/TSG.2016.2561303.

[20] Bajpai P, Chanda S, Srivastava AK. A novel metric to quantify and enable resilient distribution system using graph theory and choquet integral. IEEE Trans Smart Grid 2016;99:1. http://dx.doi.org/10.1109/TSG.2016.262381.

[21] Gao H, Chen Y, Xu Y, Liu C. Resilience-oriented critical load restoration using microgrids in distribution systems. IEEE Trans Smart Grid 2016;7:2837–48. http:// dx.doi.org/10.1109/TSG.2016.2550625.

[22] El-Sharafy M, Farag H. Back-feed power restoration using distributed constraint optimization in smart distribution grids clustered into microgrids. Appl Energy 2017;206:1102–17. http://dx.doi.org/10.1016/j.apenergy.2017.08.106.

[23] Craparo E, Karatas M, Singham D. A robust optimization approach to hybrid mi- crogrid operation using ensemble weather forecasts. Appl Energy 2017;201:135–47. http://dx.doi.org/10.1016/j.apenergy.2017.05.068.

[24] Mohan V, Singh G, Ongsakul W. An efficient two stage stochastic optimal energy and reserve management in a microgrid. Appl Energy 2015;160:28–38. http://dx. doi.org/10.1016/j.apenergy.2015.09.039.

[25] Umeozor E, Trifkovic M. Operational scheduling of microgrids via parametric programming. Appl Energy 180:672–681. doi: 10.1016/j.apenergy.2016.08.009.

[26] Khodaei A. Resiliency-oriented microgrid optimal scheduling. IEEE Trans Smart Grid 2014;5:1584–91. http://dx.doi.org/10.1109/TSG.2014.2311465.

[27] Farzin H, Fotuhi-Firuzabad M, Moeini-Aghtaie M. Enhancing power system resi- lience through hierarchical outage management in multi-microgrids. IEEE Trans Smart Grid 2016;7:2869–79. http://dx.doi.org/10.1109/TSG.2016.2558628.

[28] Gao H, Chen Y, Xu Y, Liu C. Dynamic load shedding for an islanded microgrid with limited generation resources. IET Gener, Transm Distrib 2016;10:2953–61. http:// dx.doi.org/10.1049/iet-gtd.2015.1452.

[29] Gholami A, Shekari T, Aminifar F, Shahidehpour M. Microgrid scheduling with uncertainty: the quest for resilience. IEEE Trans Smart Grid 2016;7:2849–58. http://dx.doi.org/10.1109/TSG.2016.2598802.

[30] Haddadian H, Noroozian R. Multi-microgrids approach for design and operation of future distribution networks based on novel technical indices. Appl Energy 2017;185:650–63. http://dx.doi.org/10.1016/j.apenergy.2016.10.120.

[31] Tianguang L, Qian A. Interactive energy management of networked microgrids- based active distribution system considering large-scale integration of renewable energy resources. Appl Energy 2016;163:408–22. http://dx.doi.org/10.1016/j. apenergy.2015.10.179.

[32] Nikmehr N, Najafi-Ravadanegh S, Khodaei A. Probabilistic optimal scheduling of networked microgrids considering time-based demand response programs under uncertainty. Appl Energy 2017;198:267–79. http://dx.doi.org/10.1016/j.apenergy. 2017.04.071.

Table 7 Resiliency indices of the network feeders.

With tie & 100% DG

With tie & 50% DG

With tie & 0% DG

No tie & 100% DG

No tie & 50% DG

No tie & 0% DG

F1 0.9454 0.9454 0.9454 0.6282 0.2461 0.2461 F2 0.9397 0.6821 0.1665 0.9397 0.1665 0.1665 F3 0.9468 0.9468 0.9468 0.9468 0.7207 0.2658 F4 0.6952 0.4384 0.4384 0.519 0.2622 0.2622 F5 0.7977 0.4465 0.3094 0.7977 0.4465 0.3094 F6 0.8307 0.8307 0.8307 0.7195 0.3114 0.0969 F7 0.724 0.5915 0.4204 0.724 0.5915 0.4204 F8 0.9462 0.3321 0.3321 0.8705 0.2564 0.2564 F9 0.6772 0.461 0.461 0.4214 0.2052 0.2052 F10 0.9588 0.7356 0.4309 0.9588 0.7356 0.4309 F11 0.5666 0.2029 0.1517 0.3931 0.1822 0.131 F12 0.4001 0.3169 0.2318 0.3078 0.2246 0.1395 F13 0.9357 0.9233 0.9233 0.1123 0.1123 0.1123 F14 0.8815 0.8815 0.8815 0.6268 0.4128 0.4128

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

459

[33] Mazidi M, Monsef H, Siano P. Robust day-ahead scheduling of smart distribution networks considering demand response programs. Appl Energy 2016;178:929–42. http://dx.doi.org/10.1016/j.apenergy.2016.06.016.

[34] Chen C, Wang J, Qiu F, Zhao D. Resilient distribution system by microgrids for- mation after natural disasters. IEEE Trans Smart Grid 2016;7:958–66. http://dx.doi. org/10.1109/TSG.2015.2429653.

[35] Ding T, Lin Y, Bie Z, Chen C. A resilient microgrid formation strategy for load re- storation considering master-slave distributed generators and topology re- configuration. Appl Energy 2017;199:205–16.

[36] Hosseini M, Barker K, Ramirez-Marquez E. A review of definitions and measures of system resilience. Reliab Eng Syst Saf 2016;145:47–61.

[37] Arghandeh R, Meier A, Mehrmanesh L, Mili L. On the definition of cyber-physical resilience in power systems. Renew Sustain Energy Rev 2016;58:1060–9.

[38] The White House. Presidential policy directive 21: critical infrastructure security and resilience, PPD-21, Feb. 12; 2013.

[39] Panteli M, Trakas DN, Mancarella P, Hatziargyriou ND. Boosting the power grid resilience to extreme weather events using defensive islanding. IEEE Trans Smart Grid 2016;7:2913–22. http://dx.doi.org/10.1109/TSG.2016.2535228.

[40] Panteli M, Pickering C, Wilkinson S, Dawson R, Mancarella P. Power system resi- lience to extreme weather: fragility modelling, probabilistic impact assessment, and adaptation measures. IEEE Trans Power Syst 2016;99:1. http://dx.doi.org/10. 1109/TPWRS.2016.2641463.

[41] Ma S, Chen B, Wang Z. Resilience enhancement strategy for distribution systems under extreme weather events. IEEE Trans Smart Grid 2016;99:1. http://dx.doi. org/10.1109/TSG.2016.2591885.

[42] Pinson P, Girard R. Evaluating the quality of scenarios of short-term wind power generation. Appl Energy 2012;96:12–20. http://dx.doi.org/10.1016/j.apenergy. 2011.11.004.

[43] Ray S, editor. Graph theory with algorithms and its applications: in applied science and technology. Springer: India; 2012.

[44] Fischer J, Heun V. Theoretical and practical improvements on the RMQ-problem, with applications to LCA and LCE. Combinatorial pattern matching. Berlin: Springer; 2006. p. 4009.

[45] Yuan H, Li F, Wei Y, Zhu J. Novel linearized power flow and linearized OPF models for active distribution networks with application in distribution LMP. IEEE Trans Smart Grid 2016;99:1. http://dx.doi.org/10.1109/TSG.2016.2594814.

[46] IBM ILOG CPLEX 12.6, < http://www.cplex.com > ; 2015. [47] Zhang D, Fu Z, Zhang L. An improved TS algorithm for loss-minimum re-

configuration in large-scale distribution systems. Electr Power Syst Res 2007;77:685–94.

[48] https://www.dropbox.com/s/o601t2un21olwqv/SaadatAbad%20Distribution %20Networks%20Data.xlsx?dl=0.

[49] http://www.irimo.ir/far/index.php.

S. Mousavizadeh et al. Applied Energy 211 (2018) 443–460

460