A three-stage stochastic planning model for enhancing the resilience of distribution systems with microgrid formation strategy

In recent years, severe outages caused by natural disasters such as hurricanes have high-lighted the importance of boosting the resilience level of distribution systems. However, due to the uncertain characteristics of natural disasters and loads, there exists a research gap in the selection of optimal planning strategies coupled with provisional microgrid (MG) formation. For this purpose, this study proposes a novel three-stage stochastic planning model considering the planning step and emergency response step. In the ﬁrst stage, the decisions on line hardening and Distributed Generation (DG) placement are made with the aim of maximising the distribution system resilience. Then, in the second stage, the line outage uncertainty is imposed via the given scenarios to form the provisional MGs based on a master-slave control technique. In addition, the non-anticipativity constraints are presented to guarantee that the MG formation decision is based on the line damage uncertainty. Last, with the realisation of the load demand, the cost of load shedding in each provisional MG is minimised based on a demand-side management program. The proposed method can consider the step-by-step uncertainty realisation that is near to the reality in MG formation strategy. Two standard distribution systems are utilised to validate the correctness and effectiveness of the presented model.

hardening and DGs installation can economically enhance the distribution system resilience. For this purpose, the dependency of PS and ERS should be simultaneously taken into consideration, with a resilient distribution system planning model for the line hardening, DGs installation, and MG formation.
However, resilient distribution system planning studies have mainly used the two-stage programming model to optimise investment decisions based on network reconfiguration. In these models, the first stage makes the investment decisions and the second stage minimises the load shedding costs in emergency conditions. In [6], a novel robust programming model was presented for resilient distribution system planning under worst-case scenarios. In [19], a two-stage stochastic model was presented to protect distribution systems against extreme weather events. In the first stage, the model made resilience-oriented design measures. In the second stage, the load shedding cost during the event and the damage repair cost after the event were assessed. In [20], a distributionally planning model was proposed to improve the resilience of the distribution network against two types of uncertainties. In addition, this work has not considered demand-side management to enhance the operational resilience of the system.
On the other hand, a few stochastic planning studies have been carried out on the form of the provisional MG, leading to the enhancement of distribution system resilience, while none of the planning researches considered dynamic MG formation based on master-slave control (MSC) framework [21]. In [22], a bi-level optimisation was presented to determine the boundaries for MGs and minimise the investment cost within MGs considering the imposed reliability constraints. In [23], a renewable energy planning model was presented considering the dynamic region of the MG. In [24], a two-stage stochastic optimisation framework was proposed to integrate a microgrid investment problem and operational problems in the grid-connected and islanding modes. These works did not consider the resilience concept in the distribution system planning. In [25], a mathematical optimisation model was proposed that integrates the objectives of reliability and resilience. In that model, line hardening was not utilise as an investment measure to supply critical loads during stochastic islanding. In [8], a two-stage stochastic optimisation approach was proposed to evaluate the impacts of investment decisions and various uncertainties on the system performance during and after emergency conditions. In that work, three groups of random variables such as line damage status, repair costs and load demands were considered. However, MG formation was done based on the load demand in each scenario at the second stage. It means that the load demand uncertainty directly impacts the MG formation in each scenario while it is not appropriate for the dynamic MG formation: In fact, only on the basis of the line damage, the possible MGs can be identified, and only after their formation, it is possible to take actions to supply the critical loads. In addition, that work did not consider demand-side management as an operational strategy in critical condition. In [4], a three-stage stochastic model was proposed to allocate MEGs in distribution systems. This model combined investment decisions in the PS, pre-position decision in the PRS, and real-time operation in the ERS to enhance the distribution system resilience. However, load uncertainty (LU) is important due to the uncertainty of the time of the hurricane, but that work did not consider LU for the operation of each provisional MG. In addition, line hardening and DG allocation were not considered in preventive measures. In [5], a two-stage robust programming model was formulated to improve the distribution network resilience through optimal line hardening and the formation of multiple islanded provisional MGs. In the master problem, the line hardening cost was minimised, and the subproblem involved discovering the impact of the worst N-p contingencies related to the operation of multiple provisional microgrids. According to that work, the provisional MGs are formed based on the worst contingencies, and then the optimal operation of MGs is implemented based on the uncertainty of the renewable distributed generators output. It means that line damage uncertainty (DU) for the worst-case event and the LU are not considered simultaneously for forming and operating provisional MGs. This research gap can be addressed by formulating a three-stage stochastic programming for distribution system planning, dynamic MG formation based on the master-slave concept and optimal operation of the provisional MG to improve the resilience of the distribution system. The classification of previous literature in distribution systems is summarised in Table 1, which includes resilience measures, uncertainty modelling, stochastic variables, resilience stage and operational strategies.
Thus, in the above literature review, none of the studies considered a joint scheme of resilient planning and operational strategies such as dynamic MGs formation based on the masterslave concept, network reconfiguration and demand-side management. Moreover, mathematical formulations for modelling uncertainties in MGs formation problem are the current gap in recent research studies. Therefore, in this study, a stochastic three-stage framework based on mixed integer linear programming (MILP) is introduced in order to address the shortcoming of previous studies. The key contributions and features of the proposed approach are: (1) A novel three-stage stochastic optimisation model is developed to capture the impacts of investment decisions and uncertainties on MG formation strategy. (2) The planning decisions and dynamic MG formation are integrated with uncertainties over three stages with nonanticipativity constraints, which guarantee that the MG formation decision is independent with respect to the realisations of uncertainties. (3) In comparison to the two-stage planning model, considering two uncertainties independent is near to the reality in MG formation strategy. (4) A comprehensive formulation based on a MSC technique and load control capability is proposed to improve the distribution system resilience.
The rest of the study is divided into the following sections. Section 2 describes the problem statement. Section 3 presents in detail the three-stage stochastic model for resilient distribution

CONCEPTUAL FRAMEWORK
In this section, the perceptual framework of the proposed stochastic three-stage planning model for line hardening and DG installation is presented and compared with the two-stage planning model as depicted in Figures 1(a) and (b). In the twostage planning model, the investment decisions are made in the first stage (i.e. the PS). Then, the optimal load shedding and re-dispatch of DGs in each MG are implemented in the second stage (i.e. the ERS) according to the line DU and LU. This means that both the LU and the line DU impact, in the second stage, the MG formation. However, this is far from the reality, as the formation of the MGs should be based only on the line DU because, generally, the first goal of the system operator is to re-establish the connectivity of the network in presence of damages (even accepting controlled island) and then act to supply at least the most important loads. Thus, the line damage may indicate the potential MGs, and after their formation, it is then possible to take actions to supply the main loads. Therefore, the two-stage stochastic planning model, which considers line DU and LU together, is not appropriate for implementing dynamic MG formation strategy. Moreover, investment decisions such as line hardening and DG placement cannot be effectively made with the two-stage stochastic model. As shown in Figures 1(a) and (b), two-line damage scenarios are considered based on the fragility functions of distribution components, that is, DU1 and DU2. In addition, four independent load scenarios are considered in each line damage scenario, termed LU1-1, LU1-2, LU2-1, and LU2-2. As a result, a total of four scenarios are created, represented by Moreover, the decisions in each stage pertains to the realisation of the respective uncertainty, which is indicated by the non-anticipativity constraints. In comparison to the two-stage stochastic model, the information of future uncertainty can be effectively used for the MG formation and dispatch of DGs. In the first stage, the investment decisions are defined based on all possible scenarios. Then, the dynamic MGs can be formed based on the master-slave concept in the second stage. The dynamic MG formation is dependent on the realisation of the uncertainty (DU1 and DU2) in the second stage, which is determined by the non-anticipativity constraints. In addition, the LU is considered in the operation of each MG according to the stochastic approach. For instance, if DU1 is defined in the second stage, the dynamic MG will be operated based on the realisations of LU1-1 and LU1-2. Finally, with the further realisation of LU, the DGs can be dispatched and the optimal load shedding can be employed with the operational strategies of dynamic MG formation, demand-side management and network reconfiguration to improve the resilience of distribution system against hurricanes.

MATHEMATICAL FRAMEWORK
In this section, a stochastic three-stage optimisation model for resilient distribution system planning is built, including the planning constraints in Section 3.1, dynamic MG formation and non-anticipativity constraints in Section 3.2, provisional MG operation in Section 3.3, and uncertainty modelling in Section 3.4.

Objective function and the first-stage constraints
In the first stage, the line hardening and DG installation are decided with the consideration of DU and LU. Therefore, the objective is to minimise the investment and the expected costs Stochastic formulation for resilient distribution system planning (a) two-stage model, (b) the proposed three-stage model for load shedding based on load priorities and size. It is worth mentioning that the investment schemes are shared among all scenarios. The flowchart of the proposed planning model is shown in Figure 2. The objective function is presented as follows:

FIGURE 2
The flowchart of the proposed three-stage model The binary variable Z EM m,i is employed to state whether DG m is installed at node i. Constraint (2) ensures that each node is located with a maximum of one DG. Constraint (3) states that each DG can be located to only one node. Constraint (4) indicates that the node unsuitable for DG installation should be set to zero. Constraint (5) indicates that only the DG being allocated in the PS can be installed in the candidate node.

The second-stage constraints
Generally, in the ERS, operational strategies such as network reconfiguration and dynamic MG formation are deployed considering various scenarios related to line damages. The line DUs such as DU1 and DU2 are shown in Figure 1(b). Therefore, the dynamic MGs can be formed based on each DU. In this study, based on the investment decisions in the first stage, the dynamic MGs are formed based on the MSC framework in the second stage considering the DU realisation. The MSC technique is utilised in the proposed model, which means that in each provisional MG, the frequency and voltage of the provisional MG are controlled by one DG as a master unit, while other DG units are the slave that follow the set frequency and voltage [21]. Further explanation of the master-slave concept is available in [21]. According to the above statements, topological and electrical constraints are given as follows:

Connection constraint
Due to the topological limits, each node will only belong to one of the MGs or to no MG in each scenario. Constraint (6) indicates the connection limit of each node to the MGs:

Root nodes constraints
Constraints (7) and (8) represent the status of root nodes. The node i in scenario s can be connected to MG k, only if the kth member of the set Ω B MDG is chosen as the root node:

Boundary line constraints
Constraint (9) provides the status of a line with nodes on both of its sides not belonging to the same MG in each scenario. This equation is transformed to the linear constraint based on linear methods. Hence, the boundary line constraint can be linearised using Equations (10)-(13) [5,15]: l ,k,s ≤ l ,k,s ≥

Constraints of the line and node status
Constraints (14)- (17) describe the status of a line and the damage state of nodes in each scenario. Constraint (14) connects the functional status of distribution lines with their damage status and the line hardening strategy. At first, z l,s for all distribution lines in the scenario generation phase is sampled. Then, the active status of a distribution line l,s is determined based on the first-stage decision Y l and the uncertainty realisation. Constraint (15)

Radiality constraints
In this section, the approach introduced in [21] is utilised to model the radiality of the distribution network. In this method, a fictitious network is made, where each provisional MG allows to have only one DG as a master control unit, and all other nodes are sink nodes that have loads. According to [21], the connectivity and radiality constraints can be stated by the following equations: i,r,l ,k,s = − r,k,s × r,k,s ,

Constraints of load demand
In this study, it is considered that the DSO can directly control some of the loads by utilising the demand response program. Equations (25) and (26) represent the scheduled active power for each load at the second stage considering the load control capabilities. Moreover, the load control capabilities are modelled by Equations (27)- (29):

DG constraints
The scheduled active power of DGs is restricted in Constraint (30). It can be seen that Constraint (30) is a non-linear constraint and should be transformed to the linear model. Therefore, Constraint (30) could be linearised using Constraints (31)-(34) [26]:  anticipativity constraints have been used in problems such as MEG planning, natural gas and power network expansion planning, and market clearing [4,27,28]. According to Section 2, it can be observed that i,k,s in the three-stage planning model is a decision variable, which keeps updating based on the DU. On the other hand, the MG formation decision in the two-stage planning model is represented by the variable i,k,s , which changes with the load DU. In addition, the value of i,k,s must be dependent on the realisation of DU. This study presents the non-anticipativity constraints to fill this gap. At first, the scenarios with the same DU realisation should be categorised into a particular set, defined as Ω S ′ . In this regard, it is described that DU S and DU S ′ are the DU realisations in the scenarios Ω S and Ω S ′ , respectively. Hence, if DU S = DU S ′ = DU is consented, the scenarios s and s ′ are placed into the set Ω S ′ . For example, there are four scenarios and two DU realisations in Figure 3. The scenarios S 1 and S 2 share the same DU realisation. Similarly, it can be observed that DU S 1 = DU S 2 = DU 1. Therefore, the scenarios S 1 and S 2 can be placed into Ω S ′ (DU 1). Likewise, S 3 and S 4 can be placed into the Ω S ′ (DU 2).   As depicted in Figure 3, the non-anticipativity constraint shows the S1 ∈ Ω S ′ (DU 1) and S2 ∈ Ω S ′ (DU 1) in Constraint (35). Hence, S 1 and S 2 should share the same MG formation decision. Similarly, Constraint (36) illustrates the nonanticipativity constraint that is performed for the scenarios S 3 and S 4 . r,k,s1 = r,k,s2 , ∀r, k r,k,s3 = r,k,s4 , ∀r, k -node distribution

Load shedding cost (k$)
Without DG placement With DG placement FIGURE 6 Impacts of DG placement on distribution system resilience

The third-stage constraints
With the realisation of LU in the ERS, the critical loads are restored with operational strategies of demand-side management in each provisional MG. Compared with other real-time operation models for MG formation such as in [5,8], the

Load constraints
Constraints (37) and (38) define the amount of active power supplied at the third stage. In addition, the reactive load is reduced according to the power factor as shown in Constraint (39). Also, the load controls are formulated by Equations (40)

Line flow constraints
Equations (49) and (50) denote the active and reactive power flows of the distribution lines.

Power balance in each node
Constraints (51)

Node voltage constraints
The restrictions of the node voltage magnitude and angles are modelled using Constraints (61) and (62). Constraint (63) states that the voltage magnitude of DG k as a master unit is set at the controlled value. Besides, Constraint (64) forces the voltage angle of the mentioned DG to be zero.

Load flow constraints
In this study, the method presented in [26] is used to perform power flow calculations in the distribution network. In this method, the magnitude and the angle of the node voltages are calculated based on an acceptable linear approximation. The linearised power flow constraints can be shown using Equations (65)-(67). In addition, the limits of the slack variables are defined in Equations (68)-(70) to build equality constraints that are valid when both nodes of a distribution line do not belong to the same MG:

Line damage uncertainty
In this section, a new approach, described in our up-coming study, is utilised to model the line DU. In the literature, line damage scenarios are constructed based on the intensity of hurricanes. On the other hand, due to a lack of sufficient information related to HILP events such as hurricanes, it is proper to generate line damage scenarios based on the worst-case event.
For this purpose, we consider two thresholds to generate line damage scenarios related to the worst-case event. Due to the experiences of DSOs in the face of hurricanes, they can select the acceptable thresholds for the line damage scenarios generation. Using the pole fragility function, the failure probability of the distribution line (i,j) can be computed as follows [29][30][31]: Hence, in each scenario, line (i,j) is considered out-of-service (z l ,s = 0) if the failure probability exceeds a pre-determined threshold. The procedure of line damage scenario generation is shown in Algorithm.

Load uncertainty
In this study, the Monte Carlo simulation is utilised for generating load demand scenarios. The distribution network load demand scenarios can be modelled using the normal probability density function. Then, the load scenarios with low probability are removed in a process called scenario reduction. The scenario reduction process decreases the computational time and cost of considering all possible scenarios. In this study, the forward reduction method is used to decrease the number of load scenarios [32].

CASE STUDY
The proposed MILP model is implemented to the IEEE 33node distribution system and the modified IEEE 69-node distribution system. Similar to the Texas coast [33], the occurrence of one hurricane of category 3 is considered. It is assumed that in the worst case, the duration of failures takes 15 h (between 10:00 and 12:00 AM). It is considered that the initial investment cost of DGs is 1000 $/kW and hardening cost for each pole is $6000 [34]. The length of each line is assumed to be proportional to its resistance. Therefore, we can compute the number of distribution poles in each line. Also, we consider that the span between successive poles is 45.72 m. The lifetime of resilience improvement measures is considered to be 15 years. Considering the interest rate to be 10%, the annual capital cost of purchasing and installing each measure is one-tenth of the initial capital cost. It is assumed that there are five load priorities, and the basic load shedding penalty cost is considered to be 14 $/kWh [35]. In the third-stage objective function, the parameter of load shedding cost is the product of the basic load shedding penalty cost and the load priority. Figure 4 shows the assumed multipliers of the active load profile on a typical day in the summer [36]. Also, we make an assumption that the multipliers of all individual buses are the same. The upper and lower bounds for the voltage magnitude range are taken to be 1.1 and 0.9 pu, respectively.
The line damage modelling of the 33-node test network is done under different thresholds and of the distribution line failure probabilities. First, the distribution line failure probabilities are computed based on the maximum wind speed of the category-3 hurricane and the fragility curve of the poles and lines [29,30]. Then, two thresholds of 15% and 20% are presumed, and the distribution lines with failure probabilities higher than these selected thresholds are considered out-ofservice. Based on this approach, two scenarios for the line damage status are constructed. Afterwards, 1000 load demand scenarios are initially generated using a normal distribution, and then they are reduced using the forward reduction approach to three scenarios. Therefore, six scenarios for the model are defined.
Our simulation was performed on a PC with Intel Core i7 CPU @3.20 GHz and 32 GB RAM. The proposed model was solved using the CPLEX solver under the GAMS environment with a 0.01% relative optimality gap.

IEEE 33-node distribution system
This test network has one upstream substation, six normally closed lines, and five tie lines. It is considered that controllable back-up DGs are natural gas-fired Combined Heat and Power (CHPs) with 500 kW capacities for improving the distribution network resilience. Due to the budget limitation, the total number of candidate DGs is confined to two. The nodes 8, 11, 21, 24, and 30 are selected as candidate positions for DG installation. The detailed node and line data are available in [7,37] It is worth mentioning that the realisation of DU-1 is shared among S 1 , S 2 , and S 3 . Therefore, the dynamic MG formation decision should be the same in scenarios S 1 , S 2 , and S 3 according to the non-anticipativity constraints. Based on a similar analysis, four distribution lines including 7-8, 12-13, 23-24, and 24-25 are hardened and two new DGs are allocated at nodes 21 and 24 in S 4 -S 6 with tie lines. In these scenarios, due to network reconfiguration, T1, T3, and T5 are closed, and one provisional MG is formed. On the other hand, if the tie lines are not considered in the planning model, two provisional MGs are formed to restore the loads. As can be seen from Figure 5(d), the amounts of loads restored in these scenarios are lower than the scenarios with tie lines. The planning results of the 33-node system are presented in Table 2. As shown in Table 2, without utilising the existing tie lines, the value of the objective function and total load shedding cost is more than the state with tie lines. In fact, these results indicate the importance of this flexibility in enhancing the recovery process. Moreover, we have compared the planning results when the dynamic MG formation strategy is not used in the proposed model. As illustrated in Figures 5(e) and (f), the load shedding increases without the utilisation of normally closed lines and tie lines. Consequently, considering the MG formation strategy can significantly reduce the total cost of load shedding as presented in Table 2.

Impact of line hardening
To validate the effectiveness of distribution line hardening, we investigate the provisional MG formation with various line hardening numbers from NH = 0 to NH = 4. In Table 3, the line hardening cost and total cost related to load shedding are listed. Without considering the line hardening, the line DU will cause $2,635,652.69 of load shedding cost in the distribution network. If we harden one distribution line (7)(8), the cost of loss of load will be lessened to $2,532,714.98. However, the cost of load shedding will not always reduce as we harden distribution lines. The reason is that hardened lines cannot change the formation of provisional MGs. In Table 3, when NH = 2 and 3, the line hardening schemes are the same, for example, two distribution lines 7-8 and 24-25 are hardened. These show that the cost of load shedding would be lower than that of line hardening cost. By hardening four distribution lines, the cost of load shedding will be decreased from $2,635,652.69 to $2,457,229.58, a reduction of more than 6.5%, which demonstrates the effectiveness of distribution line hardening. Obviously, the simulation results offer an effective strategy of determining the distribution line hardening by comparing the marginal benefits.

Impact of DGs
The distributed generators can supply power to connected loads through provisional MG formation when the distribution network experiences an outage and the upstream network is unavailable. To study the significance of DG in distribution systems under hurricanes, a comparison among distribution network planning with and without optimal DG placement is done when the distribution network is cut off from the upstream network. No DG units are considered in the distribution system for the state without DG. Optimal DG placement is performed by solving the three-stage model with investment decisions involving both line hardening and DG installation. As shown in Figure 6, the effectiveness of line hardening in terms of load shedding cost reduction is enhanced by DG installation, as the cost of load shedding in the optimal DG placement state is less than that of the state without DG. According to this figure, an increase in the line hardening number would not be effective when no DG units are installed in the distribution network. Additionally, with a coordinated distribution network investment solution involving DG allocation and line hardening, the distribution system is more resilient than the state without DG. Actually, when a distribution network is disconnected from the upstream network, the loads can be supplied by a DG if available. The DG units with hardened lines can form provisional MGs where the power will be supplied by the DG unit in each MG. This interesting observation demonstrates the significance of coordinated distribution system planning process of DG placement with line hardening to form provisional MGs for enhancing the resilience of distribution systems in a hurricane.

IEEE 69-node distribution system
The second case study involves a modified 12.66 kV distribution system with one feeder, 69 nodes, eight normally closed lines, and three tie lines. This is a distribution system with total active and reactive demands of 3.8 MW and 2.69 Mvar, respectively. The distribution system contains natural gas-fired CHPs as DGs with 500 kW capacity. The total number of newly installed DGs is restricted to three. The candidate positions for DG installation are 21, 30, 46, 50, and 61. The complete data about the line parameters and distribution system loads can be found in [38,39]. The optimal investment schemes resulting from solving the proposed model are depicted in Figure 7(a) considering line hardening number NH = 8 when the upstream network is unavailable. The total investment cost in the 69-node distribution system is $202,800. In Figure 7(b), there were 15 faulty distribution lines isolated by the normally closed lines in scenarios S 4 -S 6 as an illustrative example. In these scenarios, the distribution system sectionalised itself into one provisional MGs. In this provisional MG, the DG in node 11 is selected as the master unit and the other DGs are slave units. These results show that the proposed three-stage model can form dynamic MGs through MSC capability for recovering the critical loads. For the sake of comparison, the cost of load shedding under six different scenarios for the cases with and without investment schemes are provided in Table 4. The expected load shedding cost with optimal investment decisions is much lower than the state without investment decisions when the upstream network is unavailable. It can be observed that the optimal investment plans can directly lessen the economic losses during hurricanes.
For the loads with consumption higher than 150 kW, load control capabilities at five levels are assumed. In the first four levels, the load is decreased by 100 kW in each level, and in the last level, the load is shed completely. The utilisation of load control capabilities during the failure period is illustrated in Figure 8. The results of this section can be used by the DSO for the appropriate development of load control facilities. For instance, the control capabilities for the loads in 11 and 12 would not be beneficial to distribution system performance for the critical situation. On the other hand, due to the higher value of control options in the loads 61 and 64, increasing the load control facilities can be more effective from the resilience point of view.

CONCLUSION
This study presented a novel three-stage stochastic planning model for improving the distribution system resilience. The proposed method formed provisional MGs based on line damage scenarios and the load shedding cost is minimised according to LU in each provisional MG. The presented approach was implemented on IEEE 33-node and 69-node distribution systems, and the results demonstrated that the optimal line hardening and optimal DG placement can considerably improve the resilience of distribution networks under emergency situations.
Additionally, the proposed model indicated that the MG formation strategy and the topology reconfiguration can significantly reduce economic losses during hurricanes. In the future, we will focus on the resilience of combined power and natural gas distribution systems through the stochastic optimisation model. 1 if DG m is installed at node i and 0 otherwise ZP l,t,s slack variables for active power equality constraint ZQ l,t,s slack variables for reactive power equality constraint i,k,s binary variable indicating that node i belongs to microgrid k l ,s 1 if the line is active and 0 otherwise l,k,s 1 if the line in microgrid k is active and 0 otherwise i,k,t,s voltage angles of nodes P LC j,k,t,s difference between scheduled and deployed load control P DG m,k,t,s difference between scheduled and output power of DGs l ,t ,s a small amount i,k,s,m variable for linearizing (1) k,s fictional flows on the distribution line l (2) i,r,l,k,s binary variable for microgrid formation (3) r,k,s fictional supply of master DG L.s j,d,t,s scheduled block for the load controls L.dep j,d,t,s deployed block for the load controls