A microgrid formation-based restoration model for resilient distribution systems using distributed energy resources and demand response programs

118-bus distribution networks. Numerical results verified the efficacy of the proposed method as well.


Introduction
In recent years, inevitable natural disasters have become more severe and frequent due to climate change. These HR events lead to significant growth in power outage intensity and frequency. Among recent experiences of HR weather events causing power outages in distribution systems (Najafi et al., 2018;Panteli et al., 2017), hurricane Harvey can be mentioned which made landfall on Texas and Louisiana in August 2017 and left around 0.3 million of customers in a power outage throughout Texas for 14 days. The worldwide concerns on adverse impacts of HR weather events on critical infrastructures have led to the introduction of the concept of power system resilience. Resilience enhancement strategies aim at improving the power system response and restoration against HR events (Gholami et al., 2018).
Many researchers have studied distribution system restoration which aims to recover critical loads by resource rescheduling and structure reconfiguration after the incident. However during HR events, the upstream network may be unavailable and thus the distribution network should be operated in the islanded mode. In such conditions, traditional operational strategies cannot guarantee the continuity of power delivery (Hemmati et al., 2021). Distributed generators and smart grid technologies such as demand response programs can enhance the distribution system resilience by increasing the load restoration capability.
With the high penetration of DGs, MG formation is an effective operational strategy to restore critical loads as major faults occur in distribution systems. Optimal methods based on heuristic search (Bajpai et al., 2016;Sedzro et al., 2018;Xu et al., 2016;Zadsar et al., 2017) and mathematical programming (Alizadeh et al., 2020;Biswas et al., 2021;Gilani et al., 2020;Momen et al., 2020Momen et al., , 2021Sedzro et al., 2017;Zhu et al., 2020) were used for MG formation problem in the resilient distribution systems. In (Bajpai et al., 2016), a methodology based on the graph-theoretic method and Choquet integral was introduced for MG formation to quantify the distribution system resilience during extreme contingencies. In Sedzro et al. (2018), a heuristic approach was proposed that allows solving the post-event MG formation problem for

Nomenclature
Indices and Sets i,r set of buses, ranging from 1 to N bus j set of buses with control capability q set of buses without control capability sh set of buses with shifting capability (sha) cu set of buses with curtailing capability (cur) tr set of buses with transfer capability (tra) dg set of functional DGs, ranging from 1 to N dg mdg set of functional master DGs, ranging from 1 to N mdg pv set of functional photovoltaic units, ranging from 1 to N pv w set of wind turbines, ranging from 1 to N w e set of energy storage units, ranging from 1 to N e i DG set of buses with DG i MDG set of buses with master DG i sor set of buses with photovoltaic unit i wnd set of buses with wind turbine i ES set of buses with energy storage unit i From set of initial buses of line l i To set of terminal buses of line l L set of lines and tie lines, ranging from 1 to N L L nots set of lines without switch m set of formable microgrids, ranging from 1 to N mdg t index of the post-event restoration time ψ set of scenarios, ranging from 1 to N ψ θ max maximum limit of the voltage angle τ set of time slots for transferable demand response, ranging from 1 to 24  Xu et al. (2016), a resiliency-based methodology was presented to use MGs for restoring critical loads after a major disturbance. This work considers the stability of MGs and technical bounds on the transient current and voltage of DGs as well as frequency deviation as constraints of the load restoration problem. In Zadsar et al. (2017), a metaheuristics algorithm was proposed to optimally operate the smart distribution network considering MG formation in the presence of distributed energy resources. The authors neglected different demand response programs for the optimal operation of the distribution system. In Zhu et al. (2020), a MG formation model was proposed considering voltage and power loss constraints in the operation and power balance feasibility. The model considers exact power flow equations tackled by a mixed-integer linear programming (MILP). In Biswas et al. (2021), a chance-constrained optimal distribution network partitioning was proposed for identifying the best candidates for microgrid formation after HR event.
On the other hand, various types of control strategies can be deployed for DGs that may result in different operation schemes for system restoration. For example, droop-control method was widely used as a control strategy for DGs. Droop-based methods do not need communication among DGs. However, current circulation among DGs is recognized as a main drawback of droop control method. Therefore, the master-slave control concept was utilized to solve the mentioned problem. In this concept, one DG unit can control the voltage and frequency of the distribution system, which is stated as the master unit. The rest of DG units are stated as the slave units that work in current control mode. Thus, the master-slave control technique (Alizadeh et al., 2020;Gilani et al., 2020;Momen et al., 2021Momen et al., , 2020 is frequently used for MG formation problem in distribution system resilience studies. In Alizadeh et al. (2020), a bi-level optimization model based on master-slave technique (MST) was presented to boost the distribution system resilience after natural disasters considering the availability of fast-charging stations. In the lower level, the dynamic charging demand of in-service stations is determined according to the transportation network constraints and the upper level determines the boundaries of the islands for maximizing the recovered loads. In Gilani et al. (2020), an MILP model was presented to recover critical loads based on MG formation and optimal management of distributed energy resources after extreme events. In Momen et al. (2021), a two-stage optimization was proposed for enhancing the resilience of distribution systems after an HR event using distributed energy resources. A single-objective model was introduced in (Momen et al., 2020) to form MGs by means of electric vehicles and direct load control to restore critical loads. The proposed model in this work reduces the number of binary and continuous variables leading to significant enhancement in computational performance of the problem. A method was presented in Sedzro et al. (2017) to optimally form MGs aiming at recovering critical loads after a disturbance. However, these works did not consider various demand response programs including transferable, curtailable, and shiftable loads to improve distribution system resilience in extreme conditions. In addition, load restoration is just considered as a single-objective optimization in these works.
Implementing demand response programs is an effective measure for improving the resilience of distribution systems against HR events. Critical loads could be restored more efficiently as non-critical loads are shed or shifted according to contracts between customers and distribution system operator (DSO). In Shi et al. (2021), a post-event restoration model was presented to enhance distribution system resilience by considering distributed energy resources. This work considers critical loads and interruptible loads as a demand response program to enhance the resilience of distribution systems after HR events. In Bynum et al. (2019), a grid-centric stochastic programming model was proposed to employ demand response for enhancing network resilience instead of investing in transmission line hardening. However, this work considered  Res Net ψ network resilience index demand response contracts for chemical process facilities in transmission network. In Hafiz et al. (2019), a framework for resilient distribution service restoration was presented by considering the integrated control of household flexible appliances to shape the demand curve. This work did not use dynamic MG formation to restore critical loads after HR events. Also, restoration cost is not considered as an objective function for implementing demand response programs. The classification of previous resilience studies is summarized in Table 1 which includes test system, formulation, power flow (PF), solving method, uncertainty modeling, objective functions (OF), and operational strategies.
Thus, in the above discussions, none of the articles considered a joint scheme of dynamic MG formation along with different demand response programs as operational strategies in the service restoration process. In addition, there is a gap in recent studies for a mathematical formulation addressing the multi-objective optimization of load recovery and operation costs in the post-event service restoration process. Therefore, in this paper, a cost-restoration stochastic optimization is introduced as an MILP model in order to fill the gaps of previous research studies. The key contributions of the proposed approach are the following: • A novel two-objective stochastic optimization model is proposed to maximize the network load restoration and minimize operation costs considering the uncertainties of load, wind speed, and solar radiation. The proposed model examines efficient operation of MGs under real conditions after HR events.
• The ε-constraint method is used to solve the two-objective model and the fuzzy satisfying technique is deployed to choose the best possible solution. • MG formation and demand response programs including transferable, curtailable, and shiftable loads are employed for enhancing the distribution system resilience.
The remainder of this paper is organized as follows. The twoobjective optimization model is formulated in Section 2 incorporating MG formation and demand response programs. The ε-constraint method and fuzzy-based technique for compromising the optimal solution are discussed in Section 3. Numerical studies are presented in Section 4. Finally, the conclusive remarks are drawn in Section 5. Fig. 1 shows the flowchart of the post-event restoration model for distribution networks against HR events. As damaged lines and buses are determined in Stage 1, the schedulable area for microgrid formation is identified in Stage 2. Schedulable area is the set of nodes and lines not connected to the upstream network. In the presented model, the electricity distribution network is divided into several MGs using distributed generators and automatic switches to minimize the load not supplied and operation costs. Therefore, a multi-objective restoration model is proposed in Stage 3 considering MG formation and operation constraints as well as smart grid technologies. Resilience indices are then calculated to assess the distribution network performance in the restoration process. Detailed formulation of the objective functions and associated constraints of Stage 3 are presented in the following.

Topological constraints
Graph theory algorithm is an applicable technique to consider MG formation topological constraints (Ding et al., 2017;Mousavizadeh et al., 2018). In the presented model, a network graph involving all the lines and nodes is extracted. Since tie lines are considered, the network graph will contain several loops. In addition, at least one DG with the ability to control the bus voltages and frequency in each MG is required which is called the master unit. In case of multiple DGs in a MG, only one of them will be chosen as the master unit and the others will be Table 1 Survey of previous resilience studies.  (Ding et al., 2017). Thus, the topological constraints are formulated as follows.

Reference
β pos l,m,ψ,t + β neg l,m,ψ,t ≤ β l,m , ∀l, ψ, m, t In this model, each node will only belong to one MG or it will not belong to any MGs as imposed in Eq. (1). Eq. (2) models the root node constraint. Node i can be connected to MG m if the m th member of i MDG is selected as the root node. Constraints (3)- (6) indicate that if the buses on two sides of a line are not located on the same MG, the binary variable status of the line is set to zero. Eqs. (7) and (8) denote the line status and the failure state of the damaged buses. Installing remote control switches on all lines in an electrical distribution network is not economically feasible. Hence, constraint (9) indicates the switch status on the distribution lines. A spanning tree search technique (Mousavizadeh et al., 2018) is employed to eliminate unconnected loops and nodes to ensure radiality in each MG. The radiality model is formulated by constraints (10)-(16).

Load flow constraints
In this paper, the method presented in Yuan et al., 2016) is utilized to perform load flow computations in the distribution network. Using this approach, the node voltage magnitudes and angles are computed according to a linear approximation. Linearized power flow equations in  (17) and (18) represent the power balance constraints in each bus. In addition, Eqs. (19) and (20) represent the power flow constraints of distribution lines. Also, the limitations of the slack variables are defined in Eqs. (21)-(27).

Curtailable loads
Curtailable demand response allows power curtailment within a specific period while it may cause the load to rebound in the subsequent hours Song et al., 2020). The curtailable load constraints are shown in Eqs. (33)-((37). Constraint (33) indicates the net power transferred on the bus for curtailable load unit cu. Constraint (34) denotes that the sum of decreased and increased power transfers should be zero during the demand response contracts for each load unit cu. It could be implied that the increased power transfer in time t should be formerly decreased during three hours based on the predetermined coefficients. Eqs. (35) and (36)

Bus voltage constraints
Eqs. (65) and (66) show the bus voltage magnitude and angle constraints, respectively. Moreover, Eq. (67) indicates that if DG m is chosen as the master unit, the voltage of the corresponding bus is set to the desired value. In addition, the voltage angle of the corresponding bus is set to zero by Eq. (68).

Resilience indices
To assess the performance of the network and MGs, we use the network and MG resilience indices. The network resilience index could be extracted using Eq. (70), which is stated as the ratio of the recovered loads of network to sum of all connected loads to the network considering load priority weights. The MG resilience index is also calculated using Eq. (71), which is expressed as one minus the ratio of the unsupplied loads to the supplied loads in each MG, taking into account the load priority weights.

Objective functions 2.4.1. Network load restoration
In the first objective function, the performance of the distribution network against the HR event should be improved. Therefore, the first objective function helps minimize the load not supplied. In other words, the first objective function maximizes the total restored loads of the scenarios based on their priority as shown in Eq. (72):

Restoration cost
Restoration cost includes the budget required for implementing demand response programs and power supply from DGs, wind turbines, photovoltaic resources, and energy storage units as written in Eq. (73). The restoration cost should be minimized. Eqs. (74)-(81) detail different elements of restoration cost.

Modeling uncertainty of input data
Most of the models resembling real-world processes have some uncertain parameters. This section uses a scenario-based strategy to model the uncertainty of the network input data, namely the load, wind speed, and solar radiation (Farsangi et al., 2018).

Modeling load uncertainty
The normal probability distribution function is employed to model the load uncertainty as shown in Eq. (82). In this relation, variables μ and δ represent the mean and the standard deviation of the load, respectively (Farsangi et al., 2018).

Modeling the wind speed uncertainty
The Weibull probability density function is used to formulate wind speed prediction errors as shown in Eq. (83). The variables k and c represent the shape and scale indices, respectively, computed by using mean and standard deviation of wind speed (Farsangi et al., 2018).

Modeling the solar uncertainty
Beta probability distribution function, shown in Eq. (84), is used to represent the solar radiation uncertainty parameter. The variables α and β are the beta distribution parameters and are calculated according to the mean value and standard deviation of the solar radiation (Farsangi et al., 2018).
The computational burden of a scenario-based optimization model depends on the number of generated scenarios. Accordingly, it is necessary to reduce the set of main scenarios in such a way that the problem characteristics do not significantly change . The number of reduced scenarios depends on the nature of the optimization problem. This paper uses the K-means clustering algorithm to reduce the number of scenarios. More details about K-means clustering method are given in (Chévez et al., 2017).

Optimization method
The ε-constraint method is utilized to solve the proposed multiobjective optimization model. As the first step in ε-constraint method, the payoff table should be calculated. The payoff table refers to the individual optimization results of the each objective function. After the computation of the payoff table, one of the objective functions is considered as the main function while the other is considered as a constraint for the main objective function (Nojavan et al., 2018). Hence, a single-objective problem can be optimized according to constraint (85).
The single-objective optimization problem is solved for each ε (beginning with OF 2 min and ending with OF 2 max ), after which optimal solutions are derived. The set of all solutions extracted for all variations of ε are known as the Pareto optimal front of the multi-objective optimization problem. Afterward, the best possible solution is selected from among the obtained solutions by using a fuzzy satisfying technique to convert both of the conflicting objective functions to their respective normalized forms as written in (86)  .
In the above equation, μ s k refers to the membership function of the s th solution of objective function k whereas f min k and f max k represent minimum and maximum values of the kth objective function in the payoff table, respectively. Then, for each Pareto solution, the minimum membership function is identified: , ∀s Finally, the Pareto solution with the maximum value of τ s will be selected as the final solution of the multi-objective problem. From a practical perspective, the Pareto optimal front (POF) helps system operators and managers perceive the status of objective functions at each POF solution in addition to the final solution. In extreme conditions after HR events, managers may prefer to select a non-optimal solution in exchange for satisfying management issues. The decision of managers is biased by different factors including social impacts of load interruptions and economic constraint related to available budget for service restoration.

Numerical results
In this section, the proposed model is implemented on the IEEE 37node and the modified 118-node distribution systems. Simulations were performed on a PC with an Intel Core i7-700 2.4 GHz and 8 GB of memory. The proposed model was solved using CPLEX solver under the GAMS environment with a 0.01% relative optimality gap. It is worth mentioning that the proposed microgrid formation-based restoration model can be implemented for the post-event condition after any natural disasters provided that particular consequences of the event such as line outages are accurately accounted for in the model.

IEEE-37 node distribution network
This test network is configured with one upstream substation, 12 normally-closed lines, and 6 tie lines. The total active and reactive power consumptions of this network are 22.71 MW and 17.04 MVAr, respectively. Detailed information about the system node and line parameters is available in (Borghei & Ghassemi, 2021;Munikoti et al., 2021).
The required information about the location and the capacity of DGs, wind turbines, and energy storage is given in Table 2. In this test network, 2 units (DG1 to DG2) of the 3 existing DGs have the ability to act as a master unit, and hence, the maximum number of formable MGs is equal to 2. The capacity, charging and discharging rates, efficiency, and initial SOC of ES units are assumed to be 50 kWh, 50 kW, 85%, and 60%, respectively. The capacity of wind turbine and photovoltaic units are considered to be 100 and 50 kW, respectively. The associated data on load, wind, and solar power uncertainties has been extracted from Farsangi et al. (2018) and is given in Anon (2022). The characteristics of demand response contracts are presented in Table 3. The operation cost of generation and storage units and demand response programs are shown in Table 4 (Nojavan et al., 2018;Zeynali et al., 2021). Fig. 2

Table 2
Required data on generation and storage units.   Table 4 The operation cost of generation and storage units and demand response programs.  Gilani et al. shows the structure of the IEEE-37 bus network at the event onset. In this paper, it is assumed that the feeder is interrupted from the main substation. Therefore, all nodes after node 3 are in the schedulable area. To explain the effect of utilizing demand response programs and distributed energy resources on the proposed model, two cases are considered as follows: • Case 1: Single-objective operation based on minimizing load shedding, • Case 2: Multi-objective operation based on minimizing load shedding and minimizing operation cost.  In this case, the network scheduling is performed with the sole purpose of supplying the maximum load in the restoration process post the event neglecting the restoration cost. The post-event network structure is shown in Fig. 3. According to this figure, a MG has been formed under the post-event conditions. The supplied energy during the study horizon is 28,511.34 kWh and the operational cost is 2428.11 $. As shown in Fig. 3, the MG formation method has succeeded in using all three DGs in coordination via demand response programs and distributed energy resources. In this case, the network and MG resilience indices are 85.61% and 97.46%, respectively.
Given the three types of demand response schemes, Fig. 4 displays the hourly graph of load control for transferable, curtailable, and shiftable loads. As seen in this figure, a portion of the load at the peak hour is decreased by paying the corresponding cost in each of the 3 load control models. However, they intend to reduce consumption at the load peak by encouraging more consumption at non-peak hours. For example, at 17 o'clock, the shiftable, curtaible and transferable loads are 58.5, 55.9 and 170.8, respectively. Moreover, the value of load reduction in this case is 741.02 kWh. For the sake of brevity, the rest of results on generation profile of DGs and renewable units as well as energy storage scheduling are not presented for case 1.

Case 2
In this case, the scheduling is performed with the aim of minimizing load shedding and the restoration cost. In the first step, the front table is calculated for the two objective functions to compute the normal values. These front values are shown in Table 5. It should be mentioned that these values are obtained from the individual solution of each objective function.
Given the formulation of the main problem and the multi-objective optimization, the epsilon constraint method is used in this study to create the Pareto boundary. In the proposed model, load shedding is chosen as the main objective function and restoration cost is considered as the additional constraint. The number of solution points considered on the Pareto front is 11 and variation steps for the maximum level of the constraint are identical. Numerical results of are summarized in Tables 6 and Fig. 5.
According to Table 6, the minimum membership degree at Point 5 is at its maximum. The restoration cost in this state is 2058.60 $ and the load shedding is 6956.12 kWh. In addition, the network and MG resilience indices are 79.11% and 94.38%, respectively. As predicted, since the restoration cost is used as an additional objective function in this case, the resilience indices are lower than those of case 1. This limitation is near to real behavior of distribution system operators, since the budget for load restoration is not limitless. As shown in Fig. 5, the Pareto solution is the best trade-off between the two objective functions. Given the selected optimal Pareto point, the optimal network structure is shown in Fig. 6.   Fig. 7 with that of with Fig. 5 shows that fewer DRs have been used in case 2 compared to case 1, which is due to simultaneous consideration of cost and load restoration within the proposed multi-objective function optimization. In addition, the most frequently used program corresponds to the shiftable load type as shown in Fig. 7. For example, at 17 o'clock, the shiftable, curtaible and transferable loads are 41.4, 55.9 and 58.5, respectively. Fig. 8 shows the hourly graph of power generation by DGs. Accordingly at the early morning hours, the generation is reduced. This reduction is due to the low level of load in these hours while the upper bound limitation of 50% for load transferring does not allow further     utilization of DG capacity. Furthermore, all renewable units in the network are used in this case. Fig. 9 shows the hourly graph of generation by wind-based and PV units.
The charging/discharging power of storage units is shown in Fig. 10. According to this figure, storage units are charged in early hours and discharged in peak hours as rationally expected.

118-bus distribution network
The structure of 118-bus network and locations of generation and storage units are shown in Fig. 11. This network has 118 buses, 117 transmission lines, and 9 tie lines shown as dotted lines in Fig. 11. DGs at buses 17, 24, 59, 67, 76, and 107 are considered as master units with the capability of voltage-frequency control and MG formation. Moreover, 3 energy storage units have been considered for the test network. The capacity of PV, wind-based, and energy storage units are 50 kW, 100 kW, and 500 kWh, respectively. More detailed data for the network lines and buses can be found in .
The structure of the test system after the event onset is illustrated in Fig. 12. Pareto-optimal solution for the 118-bus distribution network.   Table 7. The number of solution points and variation steps for the maximum level of constraint are similar to those of the 37-bus distribution network. The results of this case are shown in Table 8 and Fig. 12.
According to Table 8, the minimum membership degree at Point 5 is at its maximum. In this case, the restoration cost, load shedding, and resilience index are 8803.498 $, 201,275.2 kWh, and 49.84%, respectively. As shown in Fig. 12, the Pareto solution is the best trade-off between the two objective functions. Given the selected optimal Pareto point, the network structure is illustrated in Fig. 13.
As depicted in Fig. 13, three MGs are formed in feeders 1 and 2 in the restoration process. In addition, feeder 3 is connected to the upstream network and buses 94, 95, and 96 are supplied from feeder 3 by closing tie line T8. The resilience index for each of the MGs 1, 2, and 3 are 0.95, 1, and 1, respectively. In addition, the network resilience index is 49.84%. The hourly graph of using load control for transferable, curtailable, and shiftable loads is shown in Fig. 14. As depicted in this figure, the most frequently used program corresponds to the shiftable load type.

Sensitivity analysis
A sensitivity analysis is conducted here to investigate the impact of significant factors on the results presented in the previous sections. The scenarios and results of the sensitivity analysis are given in Table 9. Scenario 1 is considered as the base scenario in which tie lines, DERs, and demand response programs are employed simultaneously.
According to the results of sensitivity analysis, the absence of tie lines imposes the largest increase on load shedding due to limitation for forming MGs whereas the resilience index in scenario 2 decreases around 23% in comparison with scenario 1.

Conclusion
A post-event restoration scheme was proposed for distribution systems after an HR event using MG formation, distributed energy resources, and demand response programs including transferable, curtailable, and shiftable loads. A multi-objective optimization model was proposed considering two conflicting objectives including operation cost and load shedding to improve the resilience of the electricity distribution network. Using ε-constraint method, the proposed model was solved to achieve the best Pareto front solution using a fuzzy-based method. The presented model was numerically examined on two standard test networks. According to the results, the network resilience index of 37-bus distribution system decreased around 3.2% as the restoration cost is considered in the two-objective model. In addition, due to simultaneous consideration of cost and load restoration within the proposed multi-objective function optimization, fewer DRs were used in the multi-objective model compared to single-objective one. The efficiency and applicability of the proposed restoration model was verified on 118-bus distribution system as well. The results of sensitivity analysis on this network confirmed that the absence of tie lines imposes the largest increase on the amount of load shedding. Also, as the demand response programs were used in this test system, the network resilience index increased around 12%. The simulation results showed that the model can make full benefit of smart grid facilities such as MGs formation, distributed energy resources, and demand response programs to rapidly restore the distribution system post the event. In future works,  we will focus on the resilience of energy hub systems through the stochastic optimization model after an HR event.