Resilience‐oriented service restoration modelling interdependent critical loads in distribution systems with integrated distributed generators

Funding information University of Malaya RU, Grant/Award Number: GPF049A-2018; Institut Pengurusan dan Pemantauan Penyelidikan, Universiti Malaya, Grant/Award Number: FP126-2019A Abstract The exponential increase in the frequency and intensity of high impact low probability weather-related events have pivoted the paradigm of research pertaining to power systems towards resilience. Power system is considered as a critical infrastructure, directly linked to the nation’s economy, security and health. Therefore, recent researchers have proposed several techniques to enhance the resilience of power systems. In those techniques, critical loads have been considered independent in nature and different metrics have been proposed to evaluate the resilience of the network. To enhance the resiliency, this paper incorporated distributed generators in the distribution network and critical loads are modelled interdependently. Furthermore, a novel resilience metric is proposed in this paper to evaluate the resilience of a distribution system. The proposed model is formulated as a mixed integer second-order cone programming problem and the efficacy of the proposed model is evaluated on IEEE 33and 69-bus systems. The competence of the proposed resilience metric is evaluated after comparison with existing resilience metrics.


INTRODUCTION
The extreme weather events severely influence critical infrastructures (CIs) which are directly linked to the nation's security, health and economy. Generally, CIs operate in an interconnected network. The nature of interdependency can be categorised as unidirectional, bidirectional and multi-directional. In the case of unidirectional interdependence one CI was dependent on another CI. In the case of bidirectional interdependence both CIs are dependent on each other, whereas in the case of multi-way interdependence one CI can be dependent upon more than one CIs to function properly or can have a chain-or loop-like interdependence intricate relationship. The electrical power system is also considered as a CI, and is at the epicenter of CIs; therefore it is imperative to make power systems resilient. Ninety percentage of outages in a power system originate in distribution systems (DS) [6]. Distributed generators (DGs) are one viable option for enhancing power system resilience. Resilience prowess of DGs has been proven when two microgrids (MGs) kept their connected load serviced for 2 weeks after the catastrophic damages of GEJE [7]. The aftermath of which expedited the acceleration of research into service restoration utilising DGs.
Over the last decade, extensive research has been conducted pertaining to service restoration considering DGs. DGoriented service restoration technique integrated with conventional restoration practices is proposed in [8]. Single step cold load pick up service restoration scheme utilising DGs is presented in [9]. The proficiency of DGs connected to MGs to form multi-microgrids with the objective of service restoration is investigated in [10]. Service restoration and loss reduction operation scheme for DS is proposed in [11]. In [12] a decentralised multi-agent approach for service restoration for a DS with integrated DGs is presented.
Natural calamities like hurricane Irene, Sandy and Harvey have pivoted the focus of researchers towards resilience enhancement of DS through DGs. Wang et al. [13] presented a framework for optimal operation and self-healing of DS using MGs consisting of dispatchable and non-dispatchable DGs with the objective of maximising load served. A stochastic linear programming (SLP) resilience-oriented scheduling model was introduced in [14] to mitigate the impacts of natural disasters. A resilience-oriented mathematical model for service restoration prioritising critical loads was proposed in [15]. The concept of continuous operating time was introduced to assess the availability of DG generation. A service restoration model incorporating DGs and energy management system (EMS) prioritising critical loads was introduced in [16]. Wang et al. [17] proposed a scheduling technique of networked MGs consisting of dispatchable DGs, energy storage systems (ESs) and controllable loads managed by EMS to maximise serviced load. In [18] a restoration model coordinating topology reconfiguration and MG formation considering different types of DGs to enhance resilience is introduced. A comparison between centralised and decentralised MG formation strategies for service restoration prioritising critical loads is conducted in [19].
More recent works pertaining to resilience-oriented service restoration are as follows. A sequential service restoration model for unbalanced DS under large scale power outages is presented in [20]. Multiple isolated microgrids are formed through coordination of DGs, switchable loads and switching operations. A multi-objective restoration model incorporating DGs and considering stability constraints is introduced in [21]. The objective function consists of maximising critical load restoration and minimising non-critical load shedding. The utilisation of mobile ESs in coordination with MGs to enhance DS resilience is investigated in [22]. Chen et al. [23] proposed a multi-step restoration model. The sequence of actions for ESs, DGs and controllable switches are determined to minimise unserved customers. A post-disaster restoration scheme investigating the deployment of transportable energy storage systems in coordination with MGs is presented in [24]. The potential of utilising DGs to restore a secondary distribution network in the aftermath of a natural disaster is investigated in [25]. It is evident from the literature review that the majority of researchers have proposed various service restoration models and techniques prioritising critical load advocating the resilient operational prowess of diverse sets of DGs. However, the practical interdependent nature of critical loads has yet to be investigated.
The majority of researchers have been using reliability indices such as Expected Energy Not Served (EENS), Loss of Load Expectation (LOLE) and Loss of Load Probability (LOLP) to evaluate their model's performance. These indices are not capable of encapsulating the spatio-temporal characteristics of HILP events. Presently, there is lack of standardisation of power system resilience metrics, however, there have been metrics proposed by researchers to evaluate operational resilience of power systems. Panteli et al. [26] proposed metrics to quantify the operational and infrastructural resilience of power systems. The operational resilience assesses characteristics pertaining to operational strength of the system, for example the ability to ensure the uninterrupted supply to customers or generation capacity availability in the face of a disaster. Whereas, infrastructural resilience assess physical strength of the power system, for example the ability to mitigate the portion of the system that is damaged, collapsed or in general becomes non-functional. The metrics were designed to represent each phase of the resilience trapezoid, namely, disturbance phase, post-disturbance degraded phase and restorative phase. These metrics assess the system operation in each phase accurately, however there were certain limitations. The authors did not consider variable value of restoration state, exact time stamps for system states are necessary to calculate metrics and to calculate cumulative system resilience is computationally complex. Liu et al. [27] used four metrics to quantify resilience namely, expected number of lines on outage (K), LOLP, expected demand nor served (EDNS) and difficulty level of grid recovery (G). The first three matrices are calculated individually, as they are required to calculate G. These matrices are capable of evaluating the affects of factors effecting grid recovery, but cannot calculate operational resilience. A resilience metric capable of incorporating different classes of loads is presented in [28]. The metric is a reciprocal of system loss and loss is calculated as an integral of weighted load loss in each time period. This metric is designed to capture the effect of loss of weighted (critical) load on system resilience and does not consider the system operational states to calculate resilience. A metric to assess network resilience based on system operational states was presented in [29]. This metric only focuses on system recovery and neglects system degradation. Moreover, it is not capable of considering the effect of weighted load as well.
Here, interdependency amongst critical loads has been modelled into the service restoration model of DS utilising DGs. The model emphasises on finding least loss solution prioritising restoration of interdependent critical loads. A pragmatic approach for simulation setup of the DS is adapted. DS is distributed into geographic regions, each of which is equipped with different classes of National Electric Safety Code (NESC) distribution poles. Disruptive events are setup as different category hurricanes based on the Saffir-Simpson Hurricane Scale and failure probability of poles is calculated in accordance with their fragility curves. Furthermore, a novel metric to evaluate DS operation resilience based on operation states under natural disaster is also proposed. The metric is capable of encompassing degraded and restorative phases of the system. It is also capable of incorporating the effect of post-restoration disturbance levels, which differ from original pre-disturbance levels. The main contributions of this paper are as follows: 1. Interdependent nature of critical loads along with DGs is integrated into service restoration model of DS. 2. A novel resilience metric is designed based on different operation states of DS to evaluate operational resilience.
The rest of the paper is organised as follows. Section 2 discusses the problem formulation. The restoration model is presented in Section 3. Section 4 discusses resilience metrics. Results and discussion are presented in Section 5. The conclusion is presented in Section 6.

PROBLEM FORMULATION
This section discusses the scenario-based event problem formulation of the spatio-temporal characteristics of the HILP events such as hurricanes, earthquakes etc. Problem formulation is further classified into four subsections. Section 2.1 illustrates the regional setup of the test systems, Section 2.2 describes interdependency setup between critical loads, Section 2.3 discusses the disruptive event setup and Section 2.4 illustrates the fault scenarios generated in accordance with the disruptive events.

Regional setup
A pragmatic approach to setup the event's characteristics is adopted in this study. IEEE 33 and 69-bus systems are taken into consideration as test systems in this study. In practice, a distribution system is composed of several geographic regions, therefore both test systems have been divided into four regions each, depicted by different colors. The test systems are distributed into regions based on number of buses. Regions can have overlapping distributions lines such as lines 5 and 18 in Figure 1.
In practice, distribution systems are equipped with either wooden or composite poles upon which the distribution lines are strung. Poles are divided into several classes in accordance with the ANSI O5.1 standard and NESC guidelines [30,31]. NESC classifies poles on factors such as, wind pressure, ice thickness, strength factor, load factor and grade of construction. These poles are used for three purposes: transmission lines, distribution lines and telecommunication lines. Classes 2, 3, 4 and 5 are deployed by the electrical distribution system [32]. Class 2 poles are the most robust and class 5 are the least robust based on the aforementioned factors. In this study each region of the test systems is assumed to be equipped with different classes of poles as shown in Table 1 and Figures 1 and 2, to setup different levels of susceptibility to HILP events.

Interdependency setup
Critical infrastructures (CIs) are at the center of the nation's economy, security and health. CIs are not independent rather they are greatly interdependent and mutually connected [33]. Interdependency among CIs is intricate, meaning each CI is dependent upon more than one CI to function properly and most links are reversible in nature. For example, electrical power is dependent on oil, water and telecom for fuel, cooling and communication, respectively, and oil, water and telecom are all dependent on electrical power for their functionality. Similarly, intricate interdependency exists among all CIs. Service restoration models under and after unprecedented events have extensively been researched. However, researchers FIGURE 2 Regional distribution and hurricane patterns of 69-bus system have yet to model the interdependency among CIs for service restoration. Therefore, in this study interdependency is modelled among CIs. In this study interdependency among critical loads is adopted to simulate the interdependency of CIs. Critical loads are classified into dependent and independent loads based on the relationship of their mutual functionality. For example, in Figure 1, bus 7 is modelled as a dependent load and buses 25 and 27 are modelled as independent loads. It can be assumed that the dependent load is a hospital which is dependent upon gas and sewage to function properly. Interdependent links between critical loads in test systems are shown in Figure 3 and are illustrated by blue dashed lines in Figures 1 and 2. In Figure 3a, bus 7 is dependent upon buses 25 and 22 for 60% and 20% of its functionality, respectively. Bus 7 is only capable of maintaining 10% of its functionality on its own. A similar relationship among critical loads exists in Figure 3b-d. These links are assumed as the authors did not find a precedent regarding interdependent loads.

Disruptive event setup
A distribution system should be resilient against an unprecedented event such as, floods, tsunami, hurricanes, earthquakes,

Fault scenarios in accordance with disruptive events
Distribution poles are the most vulnerable components of power systems against hurricanes. In 1989, more than 15,000 poles were destroyed by Hurricane Hugo [35]. Around 12,000 were destroyed during hurricanes Wilma and Katrina. Failure probability of distribution poles against different hazards is quantified and illustrated with the help of fragility functions and curves. Failure probability of all classes of distribution poles calculated from fragility functions and curves against wind speeds FIGURE 4 Fragility of curves of class 1 to class 7 of distribution system poles [36] is presented in [36] and is utilised to calculate the failure probability of distribution poles of all classes for both test systems in this study.
Failure probability of distribution pole classes against different categories of hurricanes is presented in Table 1 and is calculated from data based on Figure 4. The failure probability is calculated by cross-referencing wind speeds of hurricane categories with the fragility curves in Figure 4. For example, category 3 hurricane has wind speeds ranging between 111 and 129 mph, in this study the top tier speed of all hurricane categories is selected to calculate failure probability of distribution poles. Therefore, it can be seen in Table 1 that under category 3 hurricane pole class 2 has a failure probability of 10%, class 3 has a failure probability of 15%, class 4 has a failure probability of 25% and class 5 has a failure probability of 30%. As discussed earlier, regions of both test systems are equipped with the same class of poles, that is class 2 for region 1, class 3 for region 2, class 4 for region 3 and class 5 for region 4. Class 2 poles have the lowest probability of failure, that is 0%, 10% and 25% against category 2, 3 and 4 hurricanes, respectively. This reflects their high degree of robustness, whereas, class 5 poles have the highest probability of failure, that is 17%, 30% and 70% against category 2, 3 and 4 hurricanes, respectively, reflecting their low degree of robustness. Similarly, class 4 poles have a higher degree of robustness than class 5 and class 3 poles have a higher degree of robustness than class 4 poles.

RESTORATION MODEL
The following section discusses the restoration model. Following assumptions were made: (a) the distribution system under consideration is a balanced three phase AC system; (b) radial topology is maintained and (c) the operational limits voltage and current must be satisfied. The restoration model consists of power flow constraints, operational constraints, radiality constraints and supplementary constraints. The power flow constraints are illustrated by the following equations The active and reactive power balance of nodes in the distribution system are illustrated in constraints (1) and (2). Conditions of Kirchhoff 's voltage law (KVL) are represented by constraints (3)-(5). These constraints must be satisfied for every circuit of the network. Constraints (3) and (4) are adaptive representations of KVL, l j i is incorporated into them to obtain a satisfactory solution in the case when circuit ji is not operating. Conventionally, KVL does not need to be satisfied if the circuit ji is not operating; however, the restoration model requires a solution in such a case, otherwise the solution would be rendered infeasible. Constraint (6) is used to control the value of l j i and ensures the feasibility of (3) and (4). There are two binary decision variables x ji and y j in the model. The binary variable x ji represents the service status of branches; the value of x ji is 1 if the branch is in service and 0 if it is out of service. Whereas, the binary variable y j represents the connection status of nodes, the value of y j is 1 if the node is not connected meaning its active and reactive power demands are not supplied by the substations and 0 if the node is connected.
Conventionally, (5) is a quadratic equality, but is transformed into second-order conic constraint (7) by performing conic relaxation. Furthermore, (5) can be reformulated as standard second-order cone formulation depicted by (8) [37]. This transformers the non-linear problem into a second cone model in accordance with the conditions described in [38] and is valid for a linear objective function. The operational constraints are illustrated by the following equations Constraint (9) depicts the operational voltage limits of each node. Constraints (10) and (11) represent the active and reactive power operational limits of DGs, respectively. Constraint (12) illustrates the current-carrying limit of each branch. Constraints (13) and (14) limit the active and reactive power flowing in a branch within the maximum operation limit if the branch is in operation or is equal to zero if the branch is not in operation. The radiality and reconfiguration constraints are illustrated by the following equations.
Z b is the total number of buses in the system and n s is the total number of substations. Constraint (15) ensures that the total number of branches in operation is equal to |Z b | − n s , n s also dictates how many graphs the model will search for. For example if n s = 2 two graphs will exist in the solution; the first graph will contain the actual nodes connected to the grid station, whereas the second graph will contain the fictitious network of nodes not in service. Constraint (16) illustrates artificial power flow balance at each node. When y j = 0, the bus j is restored and there is no need for any artificial power flow as there is no connection of j with the fictitious node. Whereas, when y j = 1, bus j has not been restored and requires artificial power flow of value 1 to be supplied by the fictitious node. This ensures that a path must exist between j and the fictitious substation. Constraints (1), (2), (15) and (16) guarantee that the distribution system will maintain a radial topology after service restoration. Constraint (17) ensures that only fictitious substation supplies artificial flow. The artificial flow in a branch is limited by constraint (18), in accordance to its state of operation. Constraints (19) and (20) establishes that fictitious substation does not generate active and reactive power. A fictitious node (N f ) is a fictitious substation introduced in the mathematical model presented above. Fictitious node is a requirement, ensuring unnecessary switching operations do not take place to satisfy (15). The following equation illustrates the interdependency constraint.
Constraint (21) illustrates the dependency of a bus to one or more buses. In this study, the interdependent nature of critical loads has been modelled by classifying them into dependent and independent buses. For example, a hospital is considered as a critical load and it cannot operate efficiently, if its gas and sewage system (also considered as critical loads) are not operational, in other words a hospital is dependent on gas and sewage systems. In this scenario, hospital is the dependent entity and gas and sewage are independent entities. The constraint is applied on a set of dependent buses denoted by The constraint is applied separately on each dependent bus that is listed in Z d . The set Z i represents all the interdependency links between dependent and independent buses along with each link's degree of interdependence represented by q ji , that illustrates percentage of functionality between dependent and independent buses. The term degree of interdependence refers to the percentage of bus's functionality for which it is dependent upon another bus. y i represents the operational status of independent bus and y j represents the operational status of dependent bus. j represents threshold limit of the dependent bus upon which the constraint will decide, whether to shed it or keep it in service. The constraint will be forced to make this decision when an independent bus is shed meaning y i = 0, in this scenario, the constraint will shed the dependent if its percentage of functionality termed in this paper as cumulative degree of interdependence (CDI) is below the threshold value. c j is an auxiliary variable added to the constraint (21) for it to obtain a solution. When a decision is made by constraint (21) to shed the dependent load, that is y j = 1, in this scenario c j obtains a suitable value in order for constraint (21) to remain feasible. c j ensures that constraint (21) remains between 0 and j . The value of c j is controlled by constraint (22). The interdependency constraint is designed to handle various types of interdependent relations. Six such interdependency relations are illustrated in Figure 5, which will elaborate on the functionality of the constraint. It should be noted for all subsequent interdependency relations discussions in this section the threshold value is selected at 50%. Figure 5a represents a scenario in which one bus is dependent on more than one independent buses. In this scenario bus number 7 is dependent on buses 22, 25, and 10 with a degree of interdependence of 20%, 60%, and 10%, respectively. Furthermore, the circle around bus 7 indicates the percentage of functionality it can maintain on its own, that is 10 %. If bus 22 is out of service, the CDI for bus 7 is 80%, meaning that it can maintain 80% of its functionality. Moreover, 80% is above the threshold value, therefore, bus 7 will not be shed. If bus 25 is out of service, CDI is 40%, therefore; bus 7 will be shed. If bus 10 is out of service, CDI is 90% so bus 7 will not be shed.
Multi-way interdependence links are illustrated in Figure 5b. These interdependence links are similar to a chain structure. In this scenario, bus 2 is dependent on bus 27 and bus 20 is dependent on bus 5. If bus 27 is out of service, buses 2, 20 and 5 will not be shed. If bus 2 is out of service, buses 20 and 5 will be shed as their CDI is below the threshold value. If bus 20 is out of service, bus 5 will be shed. A loop-like interdependence relation is depicted in Figure 5c. In this scenario, if bus 27 is out of service, buses 2 and 20 will not be shed in accordance to CDI. If bus 2 is out of service, buses 20 and 27 both will be shed. If bus 20 is out of service, bus 27 will be shed as its CDI is less than the threshold, whereas bus 2 will remain in service as its CDI is above threshold.
Diverse one-way interdependence links are illustrated in Figures 5d,e. In Figure 5d, buses 28, 6, 16 and 19 are all dependent upon bus 14 with the same degree of interdependence. In this scenario, if bus 14 is out of service, none of the dependent buses are shed. In Figure 5e, buses 28, 6, 16 and 19 are all dependent upon bus 14 as well but with different degrees of interdependence. In this scenario, if bus 14 is out of service, buses 28 and 19 will be shed as their individual CDI is below the threshold value. In Figure 5f a bidirectional link between buses is depicted. Buses 7 and 22 are both dependent upon each other for their functionality but with a different degree of interdependence. In this scenario, if bus 22 is out of service, bus 7 will be shed. In contrast to that, if bus 7 is out of service, bus 22 will not be shed.
To summarise, the constraint can be applied on diverse set of interdependent relations, having diverse links and degrees of interdependence among critical loads. It can incorporate different values of threshold for each dependent bus. The constraint in combination with objective function will try to restore all the system load. The decision to shed a dependent critical load will only be made when its linked independent bus will be out of service. A dependent load will only be shed if its CDI is below its threshold value. The following equations illustrate supplementary constraint.
Constraint (23) guarantees that at least one branch is in service to connect the load. Constraint (24) establishes that buses connected to branches have the same operation state. Expressions (25) and (26) illustrate the binary nature of the decision variables x j i and y j . The objective function can be formulated as Subject to: (1)-(4), (6) and (9)- (26). The objective function shown in (27) consists of two terms. The first term prioritises the load to be restored. j symbolising the priority of loads can be calibrated to create various load groups. The second term symbolises power losses. Power losses are conventionally not considered as an objective for service restoration, however, it is taken under consideration in this study to force the model to select the solution with the least amount of losses. It should be noted that no individual weights have been assigned to the two terms of the objective function. The reason being that the both terms are insensitive to individual weights according to our investigations. However, inclusion of both terms is beneficial in achieving an optimal solution. If power loss is not considered as a part of the objective function, the restoration model will only search for a solution with lowest value of load shed, this usually results in a solution with high power loss. If load shed is not considered the restoration model tries to achieve a solution with zero power loss, which is only possible when all loads are shed, this renders the model infeasible. Furthermore, assigning variety of weightage in the form of percentage such that the cumulative percentage of the objective function is 100% , results in similar solution and similar value of objective function. The reason of which is that the two terms are interlinked in an intricate manner, such that power loss is constrained by the amount of load shed, as the current in the systems varies with it.

RESILIENCE METRICS
Resilience in general is defined as the ability of a system to withstand, adapt, respond and recover from an unprecedented disrupted event. There is a lack of standardised metrics for resilience pertaining to the electrical power system as discussed in [5]. However, several researchers have proposed metrics to encapsulate and assess the resilience of electrical power systems. Two existing resilience metrics have been employed and a new metric has been proposed in this study to assess the capability of microgrids in enhancing power system resilience.
Resilience metric focusing on weighted load lost under extreme weather event is proposed in [28]. The metric is capable of incorporating different classes of loads, that is critical, semicritical and non-critical. Each class of load is assigned a weight in accordance with their priority as shown in equation (30). As an illustrative example weights assigned are 10, 5 and 1 to critical, semi-critical and non-critical loads, respectively. Δp(x i ) is the total weighted load lost and P 0 is the total weighted load of the actual system before the disruptive events. Total weighted loss of the system is expressed by Equation (29), where M is the sampling number and x i is a fault scenario. Equation (28) illustrates the reciprocal relation between resilience and loss. Δp A metric to assess network resilience based on the transition states of network operation was proposed in [29]. Transition states of the network are presented in Figure 6, S o is the preevent stable operation state of network, S d is the post-event degraded operation state of the network which is a result of the disruptive event e j and S f is post-restoration stable operation state of the network. The performance of the network during these operational states is measured by the function P(t) such as P (t o ) measures the performance of network under stable state prior to disruptive event, P (t d ) measures the performance of the network in the degraded state after the occurrence of a disruptive event, and P (t f ) measures the performance of the network after the restoration process restores the system to a desired level. Network resilience denoted by R 2 is defined as the Resilience curve representing network operation states time-dependent ratio of recovery over the loss, this ratio is illustrated in Equation (31). The range of value of resilience metric R 2 at any given time t after a disruptive event is between 0 and 1. R 2 = 1 indicates that the network is fully resilient whereas, R 2 = 0 indicates that the network has no resilience.
A metric to assess operational resilience of power system was proposed in [27]. The proposed metric also based on resilience trapezoid of Figure 6. The authors proposed individual metrics for three system phases namely disturbance phase, postdisturbance phase and restoration phase. These metrics can also be combined by adding areas of two right angled triangles and one rectangle illustrated by blue lines in Figure 6 and labelled as I, II and III. Metric for phase I is representing system degradation is illustrated by equation (32). Phase II represents system after degradation; in this phase, the network operator will assess the system damage and devise a restoration plan. Phase II is illustrated by Equation (33). Phase III represents system restoration during which the network operator executes the system restoration plan. Phase III is illustrated by Equation (34). The combined resilience metric is illustrated by Equation (35).
An operational resilience metric should be capable of providing an accurate solution for different system responses as shown in Figure 7. Response 1 is an ideal case, during which system is able to fully withstand the disaster event and avoid degradation. In response 2, a system is degraded to some percentage but is not capable for restoring. For example, physical lines of some buses FIGURE 7 Types of system responses are severed and the system is not equipped with DGs, therefore, from an operational point of view they cannot be restored. In contrast to that, for response 3, a system experiences some percentage of system degradation and is capable of restoring to a level lower than its pre-disturbance level. The black dotted lines of response 4 depict a response which is similar to response 3 but has quick restoration time. Therefore, the metric should have a higher value in this case in comparison to response 3, as its response is fast. In response 5, the system experiences an initial blackout (possibly due to fault on its utility line), but the operator is capable of restoring some percentage (possibly with the help of DGs and reconfiguration). The metric's value should be zero during response 6, when the system suffers blackout (same as response 4) and is not able to restore probably because it is not equipped with DGs. Furthermore, a resilience metric should also be able to reflect systems response time.
A resilience metric focusing on the weight of loads, response time and the operational states of networks is proposed in this study. The proposed metric consists of two components, namely the recovery component and the degradation component. The recovery component, the expression left to the addition sign represents the ratio of weighted load recovered to the system's pre-event stable weighted load, multiplied with a ratio of degradation time to sum of system restoration and degradation times. The degradation component, right to the addition sign represents the ratio of the degraded system state to the system's pre-event stable weighted load.
The proposed metric (36) is also related to the performance states of the network illustrated in the resilience trapezoid in Figure 6. The metric consists of recovery component and degraded component. The recovery component assesses the percentage to which the system is restored as well as the time it takes for restoration. The degradation component assesses the percentage of system degradation. Total weighted load ΔP (t o ) representing the pre-event stable state, weighted load at degraded state ΔP (t d ) after the disruptive event and weighted load restored ΔP (t f ) after the completion of the restoration process are calculated using Equation (37). Weights of different classes of loads are represented by variables 1 , 2 and 3 . Consideration of weighted load allows the metric to capture the importance of critical loads, meaning losing a critical load will result in lower value of R 4 in comparison to that of losing a noncritical load. Moreover, the metric also includes system degradation time t d − t e and restoration time t f − t d . Ranges of both the recovery and degraded components is between 0 and 1, therefore, the expression of R 4 is normalised by obtaining the average of components. The range of R 4 is between 0 and 1, 0 indicates the least resilient system and 1 indicates highly resilient system. The proposed metric is designed to capture the degree of resilience which previous discussed metrics were not fully able to capture. The proposed metric provides accurate solution for all responses discussed in Figure 7.

RESULTS AND DISCUSSION
This section illustrates the implementation of the proposed service restoration model on test systems and discusses their results. Fault scenarios are generated in accordance with hurricanes of different categories. The simulations are carried out on a personal computer with sixth generation 2.6 GHz Intel Core i7 Quad-Core Processor and 16 GB of RAM. The proposed MIS-COP model was solved using GUROBI solver implemented on AMPL. The case study corroborates the premise that DGs can enhance power system resilience through service restoration considering interdependent critical loads.

IEEE 33-bus system
IEEE 33-bus system is taken under consideration in this study to evaluate the efficacy of the proposed service restoration method. IEEE 33-bus system is a 12.66 kV DS with 33 buses, 32 sectionalising lines and 5 tie lines. The total load demand for the system is 3.7150 MW + 2.3 MVar. The network data of the IEEE 33-bus system used in this study has been adopted from [39]. Location and size of DGs have been adopted from [40]. Three DGs have been incorporated on buses 8, 25, and 32, each of them rated at 3 MVA. All three DGs are allowed to operate in islanded mode and are assumed to be controllable. Three interdependency links among critical loads exist in the system as discussed earlier. Threshold value for every dependent load of both test system is set at 50%. The test system is distributed into four regions, each of which is equipped with a different class of distribution poles as discussed earlier. Regions 1, 2, and 3 each have eight buses each, whereas, region 4 has nine buses. All regions are equipped with both dedicated and overlapping distribution lines. Dedicated lines are geographically only situated in their assigned region, whereas overlapping lines are geographically situated between two or more regions. Dedicated lines can fail only if their respective region is experiencing external disturbance, whereas  Table 1 and number of lines per region. Based on line outages, number of line failure in each scenario are obtained through uniformly distributed random numbers. Uniform distribution is selected because the probability of failure each distribution pole, which falls under the same class is considered same. For exam-ple, failure probability of class 5 distribution poles under category 2 hurricane is 0.17, and region 4 of 33-bus test system is equipped with 11 class 5 distribution poles and each of them is assumed to have the same probability of failure.This is similar to the example of fair die in which the probability of occurrence of each value (1, 2, 3, 4, 5, 6) is same, that is 1/6. Therefore, 17% of the poles will fail in accordance to failure probability calculations shown in Table 2. Therefore, 2 lines will fail in region 4 under category 2 hurricane and line numbers 17 and 31 are obtained using uniformly distributed random numbers. Furthermore, only line faults are considered in this study, and it is assumed that either the line has collapsed physically, or faults are cleared through operation of circuit breakers either at the sending or receiving end of the line. Therefore, faulted lines are considered as open lines by the service restoration model. Two cases are compared in each scenario which is as follows 1. The system is not equipped with DGs but is reconfigurable 2. The system is equipped with DGs and is also reconfigurable In scenario 1, the disruptive event is experienced by region 4, incurring two lines faults, locations of which are shown in Table 2. Loads at buses 32, 33, and 18 having combined load of 0.36 MW are left unserved in the system after the disruptive event passes. For case 1 when there are no DGs present the system is unable to serve 32, 33 and 18 disconnected buses because all available physical connections are unavailable. Contrary to that, the system equipped with DGs was capable of 100% of the system load. Power losses for case 1 and case 2 are 103.75 and 16.21 kW, respectively.
In scenario 2, after the passing of the disruptive event, 20 buses totaling 2.055 MW of system load are left unserved. The damage inflicted on the system has also left interdependent critical loads at buses 2 and 7 unserved. The restoration model is unable to improve load serviceability employing system reconfiguration and the system remains 55% unserved. On the other hand, 97.49% of the total load was serviced using DGs in conjunction with system reconfiguration. Moreover, no critical load was lost in case 2 only 0.29 MW load was lost on buses 19 and 20. Power losses for cases 1 and 2 are 23.85 and 9.2 kW, respectively.
In scenario 3, the intensity of the disruptive event is highest incurring 8 line faults on the regions 1 and 2 that are equipped with relatively robust classes of distribution systems poles. In this scenario, the connection to the main feeder utility is severed which results in 100% blackout, hence the restoration model was unable to restore any load in case 1. The combined DG generation capacity is 79% of the total load demand, however, in case 2, DGs with help of system reconfiguration were able to restore 78.6% of the total load by forming three islands as shown in Figure 8. Green, purple and orange lines represent the electrical boundaries of the formed islands.
It can be observed that the critical load at bus 2 is unserved; this is because the generation capacity of DG at bus 25 has reached its limit. The only way to service bus 2 was to shed load on bus 27, however that is not possible as bus 27 is also an interdependent bus and is linked with bus 2. Bus 2 is dependent on bus 27 as shown in Figures 1 and 3. It is evident from this scenario that the restoration model chose to serve the independent bus and dropped the dependent bus as it was designed to do so. In this particular scenario bus 2 can be served only if DG at bus 25 is allowed to run in an overloaded state.

IEEE 69-bus system
IEEE 69-bus system is a 12.66 kV DS with 69 buses, 68 sectionalising lines and 5 tie lines. The total load demand for the system is 3.7918 MW + 2.6941 MVar. The network data of the IEEE 69-bus system used in this study has been adopted from [41]. The location and size of DGs have been adopted from [40]. Three DGs have been installed on buses 11, 49 and 61, each of them rated at 3 MVA. DGs are allowed to operate in islanded mode when the utility is not available and are assumed to be controllable. Four interdependency links among critical loads exist in the system as illustrated in Figure 2. Similar to the IEEE 33-bus system this system is also distributed into four regions, each of which is equipped with a different class of distribution poles. Regions 1, 2, 3 and 4 have 18, 17, 16 and 18 buses each, respectively. Fault scenarios are similar to those experienced by the IEEE 33-bus system. Details of fault calculation for all scenarios are illustrated in Table 2. Fault types and cases under consideration are similar to that of the IEEE 33-bus system.
In scenario 1, disruption is experienced by region 4, incurring 4 line faults, rendering 51.52% of load unserved, and location of faults is shown in Table 2. The percentage load restored is 99.97% for both cases, this is due to the fact that no tie line experienced a fault; as a result, system reconfiguration was able to restore the same amount of load for both cases. Moreover, only load at bus 20 was left unserved because all distribution lines connected to bus 20 experienced faults. However, there is a significant difference between power losses between case 1 and case 2. Power losses of case 1 and case 2 are 102.94 and 11.08 kW, respectively.
In scenario 2, disruption is experienced by regions 2 and 3 inflicting 13 lines faults, rendering 60.61% of load unserved, location of faults is listed in Table 2. 78.45% of the total system load is restored through system reconfiguration in case 1. Moreover, the critical load at bus 49 is left unserved. DGs with system reconfiguration is able to restore 98.7% of the total system load and serve all critical loads. Buses 14, 51, 49, 50, 66 and 67 are unserved in case 1 because all lines connected to them are out of service. Whereas, buses 49 and 50 served with DG at bus 49 in case 2 but the restoration model was unable to find any route to all other buses which were left unserved in case 1. Power losses for case 1 and 2 are 143.88 and 4.66 kW, respectively.
In scenario 3, category 4 hurricane is experienced by regions 1 and 2, resulting in 13 line faults, locations of which are shown in Table 2. It is evident that line 1 connecting the utility to the distribution system is also severed, as a result of which the system experiences 100% blackout. The system state cannot be improved without DGs; therefore no load is restored in case 1. Three islands are formed depicted by green, purple and orange lines in Figure 9 to serve 80% of the system load. The combined capacities of three DGs is 85% of the system load; therefore, the full potential of DGs is being utilised in this scenario. It can be observed that bus 59 is shed as DGs at buses 49 and 61 cannot restore it. Bus 10 is dependent on bus 59 for 70% for its functionality as shown in Figure 3. Therefore, when bus 59 is shed, the interdependence constraint forces the restoration model to shed bus 10. The energy that would have been used to restore bus 10 is used to restore other loads in this scenario.

Discussion
Analysis of the performance of the proposed restoration model and resilience metric is discussed in this section. The restoration model was designed with the objective of restoring maximum possible system load and choosing a solution with the lowest power loss. Furthermore, the restoration model was integrated with DGs and interdependency constrains with pre-existing network reconfiguration properties. The DGs are allowed to operate in both grid-connected and islanded mode, however, the model will only form islands when the main feeder utility is not available. The interdependency constraint is modelled by introducing a dependent and independent relationship among critical loads. The constraint is designed to first look for solution in which both critical loads are served; if such solution is not possible, then secondly if an independent bus is shed the model decides whether to restore or shed the dependent buses linked to the shed independent load. The decision is made by the interdependence constraint after investigating the CDI for each dependent bus in the system. Results demonstrate the capability of DGs to enhance the resilience of the DS. Two cases; one consisting of a DS capable of network reconfiguration and second consisting of a DS integrated with DGs also capable of network reconfiguration were taken under consideration in this study.
The system restoring prowess of the model is evident from Fig-FIGURE 9 Service restoration of 69-bus system under scenario 3 ure 10, it can be observed that the lowest percentage of load restored for case 2 for both test systems is 78.6% when 25% for the DS was under fault. It is also evident that the location, size and type of DGs plays an instrumental role in service restoration, for example in scenario 3 of the 33-bus system the restoration model had to drop bus 2, which is a dependent critical load due to the lack of generation capacity of the DG. The model did not explore the possibility of shedding bus 27 which is linked with bus 2 as it is an independent bus because the interdependence constraint is only checked if the independent bus is shed as it is designed to do so. Moreover, it can be observed in scenario 3, case 2 for 69-bus system bus 59 is shed as neither the DGs located at bus 49 and bus 61 are able to service it. Bus 59 is an independent critical bus upon which bus 10 is dependent with a degree of interdependence of 70%, therefore the interdependence constraints makes the decision to shed it as its CDI is below the threshold value. Furthermore, the DS integrated with DGs also decreases FIGURE 10 Percentage load restored for both test systems power losses in the system, as shown in scenario 1 of the 69-bus system. In this scenario both cases were able to restore 99.97% of the system load however there was a significant difference in power losses. Power loss for case 1 was 102.94 kW and 11.08 kW for case 2. Similar results were observed in scenario 1 of 33bus system power loss for case 1 was 103.75 kW and for case 2 was 16.21 kW.
There are two possibilities through which the full potential of DGs can be realised for service restoration, first the distribution system is equipped with fewer DGs with high generation capacity and second numerous DGs with less generation capacity. Each possibility has its advantages and disadvantages. The latter performs better when faults experienced by the system are geographically diverse. In such scenarios forming several MGs with small electrical boundaries can restore a higher percentage of system load as more DGs provide more routes from generation to loads. For example, in Figure 9 if two more DGs are installed at buses 27 and 34 the percentage of load restored can be increased significantly. This is not possible with a system with fewer DGs as route options are limited, however, higher capacity DG facilities possess more sophisticated control systems which result in better system operation and control which is a requirement in such extreme scenarios.
The resilience of restoration model results was quantified utilising four resilience metrics. Results of these metrics for both test systems for both cases are presented in Table 3. t d − t e is the system degradation time,which in this study is the duration of disaster event. t s − t d is the time taken for the system operator to assess damage, compute solution and start implementation. In this study it is the time taken by the restoration model to generate a solution. t f − t s is the time taken to implement the solution and achieve system restoration. In this study, this time is dependent upon number of switching operations for every case. It is assumed that meantime to operate one manual switch will take 1 h, in accordance to the assumption of [42].
The first metric investigated in this study quantifies system resilience based on the weighted load lost. However, there are three major shortcomings in this metric, first that there is no consideration of system degradation, the resilience of a system cannot be solely quantified by the load loss. Rather, resilience considers the system's ability to first withstand perturbations, then adapt itself to respond by restoring the system to the desired level. Second, the metric does not assess system degradation and restoration time. A system that restores quickly should be considered more resilient. Third, the range of this metric is very vast which makes assessing the degree of system resilience difficult. A comparison of metric calculations is presented in Table 3. It can be observed that metric values range from infinity to 1, infinity showing the system has no loss and is the most resilient and 1 indicating least resilient system with total blackout, making it difficult to assess the degree of resilience.
The second metric investigated in this study was based on the resilience trapezoid. The range of this metric lies between 0 and 1, which allows quantifying the degree of resilience unlike the previous metric. However, there are two disadvantages of employing this metric, first that the weighted load is not considered. Second that when a system degrades to a certain level and is not able to restore as depicted by response 2 in Figure 7, the metric calculation results in 0 indicating the system is not resilient at all. This should not be the case, as a system must be resilient if it is able to withstand the disruptive event without totally collapsing. On the other hand this metric deems the system fully resilient every time it is able to restore 100% of the system load as seen in scenario 1 case 2 for the IEEE-33 bus system, which is incorrect as the system should be penalised for getting degraded. Consideration of weighted loads provide a more pragmatic overview of the system performance as it reflects the loss of critical loads, therefore it should be taken into consideration. The third metric investigated in this study is also based on the resilience trapezoid of power system. The metric is capable of assessing system degradation as well as incorporating weighted loads. The metric consists of three phases each of which is calculated separately, based upon the area of that phase. Final value of the metric is the sum of individual phase results. The results of metric calculations for both test systems are presented in Table 3. There are three shortcomings with this metric. First, as the metric is based on upon calculating the area of resilience trapezoid, it results in the value of infinity if a system is degraded to some level and is not able to restore as the curve becomes infinite. The metric will render itself infinite for all system responses similar to responses 2 and 6 of Figure 7. Second, it is observed that the metric would result in a higher value if the restoration time is higher. For example, consider scenario 2 for 69-bus system. The percentages load restored for case 1 and 2 are 78.5% and 98.7%, respectively. Mean time taken to achieve restoration for case 1 and 2 are 7 and 4 h, respectively. In this scenario the metric results for case 1 should have a lower value than case 2, as less percentage of load is restored in greater time; however the results depict an opposite behaviour. Third, the range of this metric is also very vast which makes assessing the degree of resilience difficult.
The fourth metric investigated is the proposed metric in this study. The proposed metric considers weighted loads, the degraded state of the system and restoration time as well. Furthermore, it does not result in 0 if a system is not able to improve itself from the degraded state and 1 if a system is able to restore 100% from its degraded level. In fact, 0 and 1 are both ideal scenarios as a system can never be fully resilient neither can it be non-resilient. The metric is designed to normalise its value between 0 and 1. The results of both test systems corroborate the improvements. It can be observed that for scenario 1, case 2 in Table 3 the system is 3.5%degraded and 100% restored in 7 h. The value of R 2 is 1 and R 3 is 3753, whereas the value of R 4 is 0.799. Moreover, for scenario 2, case 1 the system is 55% degraded and is not able to restore in Table 3 the value of R 2 is 0 and R 3 is infinity, whereas the value of R 4 is 0.302 and for scenario 3, case 2 the system is 100% degraded and 83% restored in 12 h. The value of R 2 is 0.762 and R 3 whereas the value of R 4 is 0.244 as the system is penalised for being 100% degraded and taking 12 h to restore. Furthermore, it is evident from Table 3 that only the results of R 4 can be assimilated, results of R 1 and R 3 are ranging from infinity to 1 making it difficult to assess the degree of resilience and results of R 2 to do not depict the actual state of system resilience. The results of the proposed metric are capable of assimilating and depicting the actual state of system resilience.

CONCLUSION
The service restoration model considering interdependent nature among critical loads for DS integrated with DGs is proposed in this paper. Interdependency is modelled by establishing a dependent-independent relationship among critical loads. The proposed model is tested on IEEE 33-and 69-bus systems equipped with different classes of NESC distribution poles. Failure probability of distribution poles is calculated through the fragility curves of pole classes against wind speeds. Disruptive events are setup using the Saffir-Simpson Hurricane Scale.
Results prove the efficacy of the restoration model whilst maintaining the interdependent relationship of critical loads. Moreover, a novel metric for evaluating the operation resilience of DS was proposed in this paper. The metric was designed to overcome the shortcomings of existing metrics. The metric is composed of two components namely, recovery and degraded components, capturing the degraded and recovered operational states of the system. A comparison of results with existing metrics corroborates the metric's capability to efficiently capture the operational resilience of DS.

FUTURE WORKS
The service restoration model presented in this paper consider dispatchable DGs. However, in practice, MGs constitute of renewable DGs, therefore it should incorporate intermittent nature of renewable and non-dispatchable DGs. A time horizon-based generation profile for renewable DGs should be generated based on probabilistic approaches. In accordance to which a variable load profile should also be incorporated into the model. Furthermore, to simulate a more pragmatic approach of fault incursion, time-based fault profile should also be incorporated into the model. Reactive power demand at bus j R ji Resistance of branch (j, i) X ji Reactance of branch (j, i) I max( j,i ) Current limit of branch (j, i) S max( j,i ) Power limit of branch (j, i) V min) Lower operational voltage limit V max) Upper operational voltage limit P DG min Lower active power generation limit of DG j P DG max Upper active power generation limit of DG j Q DG min Lower reactive power generation limit of DG j Q DG max Upper reactive power generation limit of DG j n s number of substations M Big M parameter q ji Degree of interdependence for link (j,i) j Threshold value for dependent bus j Variables P ji Active power flow in branch (j,i) Q ji Reactive power flow in branch (j,i) P GR Active power generation at grid Q GR Reactive power generation at grid P DG j Active power generation of DG at bus j Q DG j Reactive power generation of DG at bus j V j Square of voltage magnitude at bus j I ji Square of current magnitude in branch (j,i) l ji Auxiliary variable to satisfy Kirchhoff's law when branch (j,i) is not in operation H ji Artificial power flow in branch (j,i) H G j Artificial power generation at bus j c j Auxiliary variable to satisfy interdependency constraint x ji Binary variable depicting the operational status of a branch (j,i). The variable is 1 if the branch is in service and 0 if the branch is out of service y j Binary variable depicting the operational status of the bus j. The variable is 1 if the bus is out of service and 0 if the bus is in service.