Prevention and Survivability for Power Distribution Resilience: A Multi-Criteria Renewables Expansion Model

Distributed energy resources are capable of enhancing grid resilience through island operations in contingency. This paper proposes a multi-year, multi-criteria generation expansion model to achieve distribution power resilience through renewable energy integration. In cases that distribution circuits or fuel lifeline are destroyed post natural disasters, wind- and solar-based generators form island microgrids to power the critical load. The goal is to determine the sizing, siting and maintenance of distributed energy resources such that system cost and power shortage in contingency are minimized. Moment methods and central limit theorem are used to characterize the spatial climate uncertainty, generation intermittency, and voltage variation. The research contributions are twofold. First, the expansion model achieves triple goals by meeting the annual load growth, reducing the carbon footprint, and enhancing the grid resilience manifested as prevention and survivability. Second, we present an improved direct zigzag algorithm to search for the Pareto solutions. Numerical examples show that the zigzag search outperforms evolutionary-based algorithms as it can identify higher-quality non-dominant solutions by exploring a wider Pareto frontier.


I. INTRODUCTION
At present physical hardening and distributed generation (DG) are the two primary approaches to enhancing the resilience of distribution grid. Adding redundant overhead lines and backup generation, fortifying poles, trimming vegetation, and preventive maintenance all fall into the hardening category [1]- [3]. The second approach for resilience attainment is to design and implement active distribution grid in which microgrid systems or distributed energy resources (DER), such as wind turbine, photovoltaic system, and battery energy storage, can operate in island mode to power the critical load post the disaster attacks [4], [5]. The global installation of DER units and microgrid capacity consistently increases, making the distribution grid more active in power supply. The study [6] shows that in 2012 alone, the world-wide new installation of DER reached 142 GW, representing 39% share of the new capacity addition of The associate editor coordinating the review of this manuscript and approving it for publication was Kai Li . that year. The rise of active distribution systems opens the way of deploying resilient smart grid via microgrid generation [7], [8].
EPRI [9] interprets and defines the grid resiliency in three aspects: 1) prevention, 2) survivability, and 3) recovery. Prevention aims to avoid or mitigate the distribution system damage. Survivability includes the use of backup, mobile or distributed generating units to maintain some basic electricity supply to critical loads. Recovery ought to provide for rapid damage assessment, prompt crew deployment to damaged sites and readily available replacement components.
This paper aims to design a resilient distribution system against disastrous weather events through the progressive integration of wind-and solar-based DER units. We propose a multi-criteria distributed generation expansion model to achieve prevention and survivability performance. In our model, prevention is achieved by allocating DER units and adopting preventive maintenance. Survivability is attained through the island microgrid operation in contingency. In case that several or all distribution lines are damaged or disconnected, renewable DER units can be clustered and form island microgrid systems to power the critical loads.
The remainder of the paper is organized as follows. Section II reviews the probabilistic DG planning models. Sections III performs the cost analysis of the DG system. Section IV characterizes the voltage variation and bidirectional power flow. In Section V a multi-criteria DG expansion model featured with prevention and survivability is formulated under generation uncertainty. In Section VI, we demonstrate the computational efficiency of the zigzag optimization algorithm in a 13-node distribution network. Section VII concludes the paper.

II. LITERATURE REVIEW
A large body of literature pertaining to DG siting and sizing is available, and readers are referred to the recent survey by Georgilakis and Hatziargyriou [10]. The main objective of DG planning is to allocate DER units in distribution networks to optimize the system performance, such as cost, reliability, power quality and carbon emissions. In these studies, grid resilience features such as prevention and survivability are not explicitly incorporated in the DER siting and sizing model. Recently, grid resilience is emerging as a key criterion to measure how a power system can survive and withstand under extreme weather or unusual events. Panteli et al. [11] use a set of severity risk index to determine the application of defensive islanding of power grid in extreme weather. It boosts the grid resilience by splitting the network into stable and self-adequate islands in order to isolate the components with higher failure probability, especially cascading failures. Kwasinski et al. [12] compare the availability of different DER units upon natural disaster onslaught. Their study shows that wind turbine (WT) and solar photovoltaic (PV) outperform fuel-based generators because the lifeline of fuel supply is likely damaged after the harsh weather or earthquake. Bie et al. [13] propose a quantitative framework along with several metrics (i.e. LOLP and EDNS) to evaluate power grid resilience, extending the scope of reliability measures to resilience measures. Ma et al. [14] propose a resilience-oriented design framework via line hardening, backup DG, and use of automatic switches. A two-stage mixed-integer stochastic optimization incorporating spatiotemporal correlation uncertainties is formulated to minimize the cost of investment, loss of load, DG operation, and damage repairs. Liu et al. [15] use Monte Carlo simulation to compare the grid resilience between network hardening, topology reconfiguration and microgrid operations based on IEEE 30-bus and 118-bus systems. Their study shows that topology reconfiguration and microgrid operations could outperform network hardening under extreme events.
To continuously supply electric power under extreme weather conditions, we propose a distributed generation expansion model by minimizing the maximum nodal power shortage when all or certain distribution lines are disconnected post disaster attack. The contributions of this paper are highlighted in following: VOLUME 8, 2020 • The proposed DG expansion model achieves triple goals: meeting the annual load growth, reducing carbon footprint, and enhancing grid resilience manifested as prevention and survivability. The novelty of our work is to incorporate prevention and survivability features into the traditional DG expansion process.
• Island microgrid operation is designed into the distribution system through onsite wind and solar generation. Hence, it can avoid cascading failures resulted from common cause failures, such as damaged fuel supply infrastructure.
• We present an improved direct zigzag search algorithm to solve the multi-criteria optimization model by exploring a wider Pareto front.

III. SYSTEM RESILIENCE AND COST ANALYSIS
A. THREE LEVELS OF RESILIENCE Figure 1 depicts five resilience states that evolve with time upon the extreme event attack: reliability, degradation, stabilization, recovery or repair, and renewed state. Let P(t) be the power at time t, three different levels of grid resilient behavior are explained below. Level 1: If WT, PV and battery are adopted as DER in the distribution grid, the grid is able to meet the load with no significant performance degradation provided the battery capacity is large enough. Such design is referred to as 'fully' resilient system as shown by path AD.
Level 2: If only WT and PV are adopted, P(t) is likely to decline due to their intermittent generation and the potential loss of substation. Path AB represents the degradation process. The design objective is to ensure that the average power from WT and PV remains at P d . In that way, the critical loads within a node still can be energized from time t 2 , At t 3 the system starts to recover and eventually return to P 0 at t 4 . Since critical loads are energized in [t 1 , t 4 ], the grid operating with path ABCD is also treated as a resilient system. Level 3: If there is no WT and PV or other DER units, P(t) will continue to drop after t 2 . The grid power eventually crashes at t 3 if all the lines are disconnected or the substation is damaged. Therefore, the system is considered as a non-resilient design.
Our model is presented at Level 2 where the distribution grid may lose certain amounts of power post the extreme event. Since WT and PV units could accelerate the power restoration by reducing the recovery time, this is another value-added resilience feature of onsite generation. The restoration becomes more difficult and takes a longer time if the grid is completely down.

B. SYSTEM COST ANALYSIS
The cost of a DG system comprises of four main items: installation, operating and maintenance, carbon credits, and equipment subsidies. Let i be the index of available DER types for i = 1, 2, . . . , m. Further let x ijt ∈ {0, 1} be the binary decision variable representing DER type i is installed at node j for j = 1, 2, . . . , n in year t for t = 1, 2, . . . , T . The present worth of DG cost over a period of T years, C eq , can be estimated by where a ij is the capacity cost ($/MW) of DER type i at node j, and P c i is the power capacity for DER type i. Parameter r is the discount factor. Though wind and solar generation does not require fuels, the costs of leasing the land for equipment placement is treated as part of the capacity cost.
WT and PV are prone to field failures. System reliability and availability can be improved through preventative maintenance (PM) by pro-actively inspecting and replacing aging components. Age-based PM is perhaps the most widely used maintenance policy in power industry because of the scheduling flexibility [16], [17]. The cost rate in [17] is adopted to estimate the maintenance expense. Let C mn be the maintenance cost of the DG system over T years, then where τ i is the decision variable representing the PM interval, and u i (τ i ) is the maintenance cost rate of DER type i during a PM cycle. Here t a i is the annual operating hours of DER type i, and z ijt = t k=1 x ijk is the cumulative DER units of type i installed at node j through year t. Note that F i (y) and f i (y) are, respectively, the Weibull cumulative distribution function and the probability density function, and η i and β i are the scale and shape parameters, respectively. R i (τ i ) is the reliability of DER type i, and R i (τ i ) = 1-F i (τ i ). Weibull distribution is widely used to model the failure probability of WT and PV units [18], [19]. Finally, c f i is the cost of a corrective maintenance action, and c p i is the cost of a scheduled maintenance action on generation type i.
The operating cost of the distribution grid is the expenses associated with the routine equipment checking, line inspection and trimming of vegetation. Based on [20], The operating cost, denoted as C op , is assumed to be proportional to the energy generated where b ij is the operating cost per MWh of electricity generated from DER type i at node j. The justification of Equation (4) is that tear and wear on turbines, solar PV, and distribution circuits likely increase with energy generation [20]. Some studies treat annual operating and maintenance cost as a fixed quantity, but others prefer to use cost per MWh (e.g. $/MWh). In this paper C op is assumed to be proportional to the energy delivered. Carbon credits and equipment subsidies aim to incentivize the adoption of WT or PV units. Let c ij be the credits given to the DER type i at node j whenever one MWh of renewable energy is produced. Also let d ij be the subsidies per MW capacity of DER type i at node j. The present worth of total incentives received by the DG system over T years, denoted as C in , can be estimated as The incentives can be treated as the cost savings for renewable DG system, hence c ij <0 and d ij <0 in equation (5).

IV. VOLTAGE VARIATION ANALYSIS A. BIDIRECTIONAL POWER FLOW
Though quite a few direct current (DC) transmission lines are now in place worldwide, the majority of distribution grids still operate in alternating current (AC) mode. To maintain the voltage stability, capacitor banks are required and used as the Volt/Var compensators in AC grid. Since electric vehicles, battery storage, and solar PV use or generate DC power, recently the design of DC distribution systems is receiving much attention in academia and industry. A major advantage of using DC grid is its simplicity in computation and control. Comparative studies have been performed to justify and validate the use of DC power flow model [21], [22], and their conclusion is that DC power flow model is reasonably accurate to represent AC power flow systems in strategic planning, such as long-term generation expansion.
We adopt the DC power flow model from [23] to demonstrate the multi-criteria DG expansion model on the basis that the planning is performed at the strategic level for multiple years. By referring to a single distribution line in Figure 2(a), the DC current in link j is where V DG is the nominal voltage of the DG system, P k and D k are the DER power and the load of node k, and J is the total number of nodes along the line. For instance, J = 4 in Figure 2. For a single distribution line, the index j for link is identical to the upstream node index. Similarly, the nodal voltage is given as with Note R l is the link resistance between nodes k and k + 1. Both equations (6) and (7) allow for estimating the line current and nodal voltage based on D k and P k . For a DG system, the unidirectional power flow condition always holds as long as the DER power is absorbed by the local node. The power flows reversely only if the sum of DER power exceeds the total demands of the downstream nodes. Both equations (6) and (7) are applicable to unidirectional and bi-directional power flow.

B. VARIATION OF VOLTAGE AND CURRENT
Note that I j and V j in equations (6) and (7) are random variables due to the uncertainty of P j and D j . Since the distribution function of I j and V j are difficult to derive, we propose a moment approach to characterizing their stochastic behavior. Let E[I j ] and Var(I j ) be the mean and variance of the current of link j connecting nodes j and j + 1. We have where, µ D k and σ 2 D k are the mean and variance of D k of node k. Similarly, µ P k and σ 2 P k are the mean and variance of P k of node k. The mean and variance of the voltage at node j, denoted as E[V j ] and Var(V j ), can be obtained as In next section, we show that equations (9)-(12) are the foundations in controlling the thermal stress and voltage variations under generation and load uncertainty.

V. MULTI-CRITERIA DG EXPANSION MODEL A. STOCHASTIC EXPANSION MODEL
In this section we formulate a multi-year, multi-criteria DG expansion model based on the voltage and current analysis in equations (9)- (12). The model captures spatial and temporal variations, including nodal voltage, power intermittency, load growth, and environmental concern. Subject to supply reliability, power quality, thermal stress, and emission constraints, the stochastic multi-objective optimization problem (MOOP) jointly minimizes the system cost and maximize the power resilience. Let x represents the DER sizing and siting of the nodes; and τ stands for the maintenance time of DER units. The MOOP, denoted as Model P1, is formulated as follows Model P1: Min: Subject to: Pr{P t (x) < D t } ≤ α 1 , ∀t.
In Model P1 the objective functions and the constraints involve intermittent power P ijt (ω) and uncertain load D jt (ξ ), where ω represents random wind and solar irradiance, and ξ is the demand uncertainty. Objective function (13) minimizes the expected present worth of the DG system. Objective function (14) maximizes the resilience of the grid by minimizing the maximum nodal power shortage if adjacent links are destroyed in extreme event. This ensures that the grid has the ability to withstand extraordinary and high-impact low-probability events. A typical example is the N-90 failure in Hurricane Sandy as opposed to the classical N-1 reliability design criterion [24].
There are seven constraints in Model P1. Constraints (15) is the reliability criterion defined as the loss-of-load probability (LOLP). Constraint (16) defines the scope of the voltage variation, where V max and V min are the upper and lower voltage limit, respectively. Constraint (17) (20) is the cumulative installation of DER type i at node j from year 1 to t. Constraints (21) simply states that τ i is nonnegative, and up to one DER unit of type i can be placed at node j in year k. This constraint, however, can be relaxed if more than one DER unit can be installed at the same node each year.
The resilience in Model P1 is manifested as prevention and survivability. To protect the grid operation against extreme weather, we rely on DER integration and preventive maintenance strategy. This is reflected by constrain (18) where the green energy coefficient λ t increases with DER installation. The survivability is explicitly incorporated in objective function (14) by minimizing the maximum nodal power shortage in contingency.

B. DETERMINISTIC COUNTERPART
It would be difficult, if not impossible, to directly search for the optimal x and τ given the random behavior of P ijt (ω) and D jt (ξ ). To make P1 more tractable, we convert it into a deterministic counterpart as follows.
Model P2: Min: g 2 (x, τ ) Subject to: x ijk , ∀i, ∀j and ∀t (29) x ijt ∈ {0, 1}; τ i ≥ 0. ∀i, ∀j, and ∀t Both P2 and P1 have the same objective functions, but constraints (24)- (27) are now in deterministic forms converted from corresponding chance constraints (15)- (18). The detailed conversion is provided in Appendix. The transformations are made based on the normality assumption on system cost, voltage variation, and thermal stresses that are commonly used in existing works [25], [26]. Note that Z-score in constraint (24) represents the Z-value of the standard normal distribution. For instance, in constraint (24) Z 0.95 = 1.64 given α 1 = 0.05. Similar interpretations can be applied to other constraints. In our case, P t (x) is the aggregate DER power of n nodes. According to the central limit theorem, if nodal power and load have limited mean and variance, P t (x) and D t tend to be normally distributed if n is large enough.

C. SOLUTION ALGORITHM
We modified a newly developed multi-objective optimization algorithm, called direct zigzag (DZZ) algorithm [27], and use it to find the Pareto solution for P2. The results are compared with a popular GA-based MOOP method, i.e. Non-dominated Sorting Genetic Algorithm or NSGA-II [28].
The key idea of the DZZ method is to search around the Pareto front efficiently by employing local neighborhood search procedures. The efficiency of a zigzag search is due to the expectation that an existing solution path will effectively pass through the Pareto front. But how to identify such a tight solution path close to the Pareto front is crucial for the performance of DZZ. Figure 3(a) illustrates that a DZZ search is performed on a deterministic MOOP minimization problem (i.e. DMOO) with two objective functions f 1 and f 2 . DZZ first employs a direct search (derivative-free optimization based on Hooke Jeeves and neighborhood search) to find the first Pareto optimum x 1 with the minimal f 1 . DZZ then iteratively performs zig and zag search to find a sequence of Pareto solutions x 1 , x 2 , . . . , x 9 . As shown in Figure 3(b), for a bi-objective optimization with two discrete variables, a zig search is simply a neighborhood search around the current Pareto solution, say x, and returns a new and feasible solution x 0 . A zag step subsequently searches neighboring points by finding x 1 0 and x 2 0 around the current candidate point for non-dominated solutions with smaller f 2 values. As shown in Figure 3(a), a zag search stops when no smaller f 2 can be found in the neighborhood of the current candidate solution.
Since Model P2 belongs to the mixed integer programing problem, we tailored DZZ to account for mixed type of decision variables x and τ . Note that the decision variables of installing DER units at nodes are binary, the pattern search in the Hooke Jeeves direct search is not applicable; nevertheless, pattern search for the continuous variables τ is still efficient. Thus, we apply pattern searches after every Hooke Jeeves exploration step for continuous variables.
Another enhancement of DZZ in this work is to use a broader local neighborhood search for Hooke Jeeves exploration step. The traditional Hooke Jeeves exploration searches the solutions in the smaller neighborhood N 1 (x) = {x : |x − x |V ≤ 1} around the current candidate solution x. To find the solutions along the ''valley'' that is not necessarily aligned with coordinates, we consider evaluating a larger neighborhood N 2 (x) = {x : |x − x |V ≤ 2}. When the decision variables are binary, each solution x ⊂ N 2 has exactly one or two elements (coordinates) with 0 or 1. For example, if x = (0, 1, 1, 0, 0, 1, 0) and x = (0, 0, 1, 0, 1, 1, 0) ∈ N 1 (x) or x = (0, 1, 1, 1, 0, 1, 0) ∈ N 2 (x). Hence the modified DZZ searches a much larger neighborhood for each of its Hooke Jeeves exploration step.

VI. NUMERICAL EXPERIMENTS A. NETWORK TOPOLOGY AND LOAD CONDITION
The topology of the radial distribution network is given in Figure 6 with n = 13 nodes interconnected by 12 links. This network is an expanded version of the 9-node system originally in [29], [30]. Note that prior to DER integration, a 30 MW substation is installed at node 1 with the nominal voltage V DG = 33 kV.
We plan to expand a DG system for a period of three years (i.e. T = 3). The demand of 13 nodes is characterized by the mean and standard deviation, yet no assumptions are made  regarding the demand distribution. Table 1 lists the mean and standard deviation of the nodal demands over three years. The aggregate mean and variance of the distribution system are also obtained and listed in the last row.
The load profile of 12 links can be estimated based on the demands of its downstream nodes. Table 2 shows the mean and standard deviation of the load of twelve links in Year 1. Similar results can be obtained for Years 2 and 3. In the table, R l is the resistance of the link between two adjacent nodes, and I max l is the maximum current that link l can accommodate. Since R l and I max l represent the physical properties of the distribution lines, they usually do not change over the planning horizon.

B. PARAMETERS FOR WT AND PV SYSTEMS
Assume five types of DER units are available for the placement in the 13-node distribution network. These are 1.5 MW WT, 2 MW WT, 2.5 MW WT, 1 MW PV, and 1.5 MW PV. Data associated with equipment reliability, maintenance, and repair are listed in Table 3. These values are obtained based on the WT and PV reliability analysis reports in [18] and [31].   DER unit is assumed to be h = 20 years, and the loan interest rate is r = 5% compounded annually. Table 4 summarizes the cost items associated with equipment, operation, carbon credits and subsidies for DER. In this example, only PV receives carbon credits to compensate the high installation cost, and no equipment subsidies are given to WT and PV. All cost items are node-independent, meaning a ij = a ik for given i and j = k, so do b ij and c ij . The unit of a i and d i is $/MW, and the unit of b i and c i is $/MWh. The last three columns define the characteristic wind speeds of WT power curve. Note that v c is the cut-in speed, v r is the rated speed, and v s is the cut-off speed. The values of these critical speed are based on the report from NREL [32].
In Table 5, the nodal wind speed follows the Weibull distribution, and c w and d w are the corresponding scale and shape parameter, respectively. Solar irradiance at each node follows the beta distribution where a s , b s and s m are the distribution parameters. Given the wind and solar distributions, methods for computing the mean and the variance of WT and PV power output are available in [25], [33]. The substation has zero power variance because the electricity is supplied by centralized power plants. To maximize the annual energy throughput, all PV panels are oriented towards the south if in the northern hemisphere. The PV efficiency is assumed to be η s = 0.20, and the average operating temperature is T o = 45 • C. Unlike WT that operates in 8,760 hours a year, the annual operating hours for PV is 4,380 hours on average. The rated power of PV module is 120 W/m 2 under clear sky with direct sunshine.
Other parameters associated with the planning model are specified as follows: α 1 = 0.001, α 2 = 0.9, α 3 = 0.9, and α 4 = 0.9. In addition, we set V min = 0.95 V DG and V max = 1.05 V DG , meaning nodal voltage varies within ±5% of its nominal value. The GEC target is λ 1 = 0.1, λ 2 = 0.2, and λ 3 = 0.3 over three years, respectively. To mimic the elimination of conventional generation, we sequentially reduce the substation capacity from 30 MW to 25 MW in t = 2 years, and to 20 MW in t = 3 years.
We solve the example problem using DZZ and NSGA-II algorithm, respectively. We choose NSGA-II as the benchmark out of two reasons. First, there are only very few methods suitable for solving discrete MOOP problems, and NSGA-II is generally available and applicable for MOOP problems with binary variables and stochastic constraints. In fact NSGA-II is often used to solve MOOP problems arising from many different fields including power engineering [34]- [36]. Second, we do not consider other MOOP methods because, as studied in [27] and [37], a zigzag search algorithm has been compared to a few other MOOP methods including pareto front approximation [38] and the normal-boundary intersection method [39]. The results in [37] show that the zigzag search is more efficient than some recently developed methods, including MO-COMPASS algorithm, standing for multi-objective convergent optimization via most-promising-area stochastic search [40]. Though there are other bio-inspired evolutionary algorithms, such as particle swarm optimization [41], they are mostly developed for continuous problems and typically slow in convergence.
The total number of design variables for the 3-year planning model is 200, including 195 binary variables and 5 continuous variables. DZZ algorithm and NSGA-II method are both implemented in Matlab computing environment. The numerical experiments are carried out on a PC workstation with 8G RAM and a 3.4GH CPU. We investigate two cases. In the first case, all nodes are assumed to be fully disconnected under catastrophic weather event. In the second case, only certain distribution lines are damaged post the extreme weather. We run four independent replications of each algorithm on the 13-node distribution network in Figure 6. The optimization results are compared and analyzed in detail below.

C. RESULTS FOR A FULLY DISCONNECTED NETWORK
In this case, it is assumed that all distribution lines are damaged post disastrous weather event. Figure 4 shows the NSGA-II optimization result for this case over a 3-year horizon of renewables integration. The Pareto front of 24 non-dominated solutions appears to be aligning along a smooth curve. The population size is 500 and the number of generation is 50 for NSGA-II.    Starting from this solution, DZZ applies a series of zigzag search to identify a solution path that is closely along the Pareto front. In that way DZZ can avoid evaluating many feasible solutions that are far from the Pareto front region. The optimal solutions cover a wide range of cost from $62M to $230M, larger than the NSGA-II frontier. Both figures show the maximum power shortage post extreme events decreases from 10 MW to zero with the increasing system investment. Table 6 summarizes the performance comparisons between DZZ and NSGA-II.  Figure 6 shows the detailed power generation and demand at every node for a particular Pareto solution found by DZZ. Optimal sizing, siting and maintenance of WT and solar PV VOLUME 8, 2020  are listed in Table 7 with the bi-objective values (g 1 , g 2 ) = ($125M, 3.88 MW). It shows that some nodes can produce more power than their demand, and the surplus power indeed flows to their neighboring nodes, either downstream or upstream. For instance, the DER units at Node 10 generate 4.1244 MW while its demand is only 2.1840 MW. Therefore, an amount of power 1.9404 MW flows along the upstream power line (i.e. reversely) to Node 9. Node 3 is another case of which the surplus power goes to the downstream Node 4 and to the upstream Node 2. For Node 3, the power line between Nodes 2 and 4 has bidirectional power flow.
The power resilience in Figure 6 is manifested as prevention and the survivability which are reflected in the allocated DER power at each node. Take example of Table 7, the optimal allocation of WT and PV capacity to each node is listed along with the equipment maintenance interval. The nodal mean power shown in Figure 6 actually demonstrates the degree of survivability of the network. For instance, the mean power of Node 2 is 3.2882 MW. If Node 2 is isolated because its adjacent distribution lines are damaged, the local power shortage is only 3.644-3.2882=0.3558 MW.
A main advantage of the non-dominant solution set over a single objective optimization lies in the fact that the decision-maker can choose a best comprised solution by trading the cost with the network survivability. In Model P2, the capability of network survivability is manifested by the Pareto optimality front in Figure 5. If the network planner is willing to invest more budget on distributed power capacity, the maximum nodal power shortage can be further reduced. Figure 7 compares the optimization performance of NSGA-II and DZZ algorithms applied for this example case with multiple replications. The plotted solutions (in the objective function space) from both algorithms are non-dominated solutions among all the evaluated solutions. For comparison purpose, we run both algorithms four times independently with different initial solutions. For all the four runs of DZZ, 2,500 to 3,000 solutions are evaluated, while NSGA-II evaluates about 25,000 solutions for each of its four runs. From the results, we see that all four DZZ runs outperform NSGA-II by locating better trade-off solutions and larger coverage of the Pareto front.
Based on our numerical experiments, the total computing times for both algorithms are approximately proportional to the number of feasible solutions evaluated. The computing time in addition to function evaluations in DZZ algorithm is negligible because it employs neighborhood search. In other words, both DZZ and NSGA-II spend the computing time mainly in evaluating the objective functions. Since a much small amount of feasible solutions are evaluated by DZZ compared to NSGA-II, the DZZ requires less amount of computational time.

D. PV COST AND CONFIDENCE LEVEL
Note that the solution listed in Table 7 indicates that PV is not favorable compared to WT. This is because the capacity cost of PV is much higher than those of WT. Suppose that PV cost is reduced to $1.5M/MW and $1.0M/MW, respectively, for 1 MW and 1.5 MW PV systems, optimal solutions show that more PV units will be chosen and installed. Table 8 shows a Pareto solution obtained under the reduced PV costs. There are 14 PV1 and 17 PV2 units are now installed as opposed to only 1 PV1 and 1 PV2 in Table 7.
Since the confidence levels in Equations (15)-(18) also influence the system cost and the nodal power shortage. By referring to Model P1, α 1 represents the LOLP, and the benchmark value in our numerical study is α 1 = 0.001.  A smaller α 1 is always preferred as it results in a lower probability of power shortage. By looking at Equation (15), reducing α 1 requires more installation of WT and PV units locally, hence the system cost g 1 (x, τ ) will increase. We solve Model P2 by reducing α 1 to 0.00075 and 0.0005, respectively, whereas other conditions remain unchanged. The corresponding non-dominant solution sets are presented in Figure 8 along with the Pareto front of the benchmark solution. Other confidence levels like α 2 , α 3 and α 4 , represent the nodal voltage, thermal stress, and renewables penetration, respectively. If one increases any of these confidence value, it is anticipated similar pareto frontiers like the case of α 1 will be obtained.

E. RESULTS OF PARTIALLY DISCONNECTED LINES
Extreme events breaking down all the distribution lines are among the worse-case scenarios, but the probability of such occurrence is very rare in reality. We consider several realistic disruptive events due to severe weather or other natural disasters that may cause some lines or components to fail. To account for varying levels of potential disruptions, for the studied 13-node network in Figure 6, we choose three representative damaging events corresponding to catastrophic, disruptive, and low disruptive scenarios respectively. More specifically, we consider three arbitrary scenarios under different levels of disruptions: s 1 = {(1, 2), (3,4), (1,5), (5,6), (1,9), (1, 11)} has six links broken with probability (1,5), (6,7), (1,11)} has four links broken with probability p 2 = 0.3; and s 3 = {(1, 5), (1,9)} has two links broken with probability p 3 = 0.6. Note a nodal pair in s k for k = 1, 2, and 3 represents a particular link. Overall, we want to minimize the expected total power shortage for all the network nodes. In this setting, the second objective function g 2 in Model P2 becomes: where S ={s 1 , s 2 , s 3 } represents the set of three disruptive scenarios being considered, and each scenario represents a damaged network of disconnected clusters of nodes. We run Model P2 on the example problem with the consideration of three disruption scenarios. Figure 9 shows the solutions evaluated using both optimization algorithms. Note that NSGA-II needs to evaluate much more solutions (∼20,000) that are all inferior to those evaluated (∼2,000-3,000) in DZZ. The results for this case also indicate that DZZ outperforms NSGA-II again in terms of computational speed and solution quality.

VII. CONCLUSION
This paper proposes a multi-year, multi-criteria distributed generation expansion model to minimize the cost and maximize the power resilience against disastrous weather events. A modified direct zigzag algorithm is adopted to search for the optimal sizing, siting, and maintenance of DER units. We enhance the distribution resilience in terms of grid prevention and generation survivability in contingency. The prevention is achieved by increasing the distributed generation power, reducing the substation capacity, and implementing preventive maintenance program. The survivability is achieved by minimizing the maximum nodal power shortage across the entire distribution system. The proposed planning model is applicable to distribution networks where lines are either fully or partially destroyed. In either case, wind-and solar-based DER units are split into island microgrid systems to power the critical loads. Since wind and solar generation does not rely on fossil fuels, this is advantageous over fuel-based generating units as the fuel lifeline is likely to be damaged in the aftermath of catastrophic events. Hence the risk of cascading failures could be largely mitigated or avoided using onsite renewable energy. The variations of load, voltage and thermal stresses are characterized by chance constraints through the central limit theorem.
Two observations are obtained from the numerical example. First, the direct zigzag search is able to explore a wider frontier and returns a set of non-dominant solutions with higher quality compared with NSGA-II. Second, the testing network indicates that 30 percent of renewable penetration can be achieved sequentially in three years by concurrently injecting WT and PV units, and reducing the substation capacity. There are several potential extensions of this work. For instance, the model can be formulated and solved in a large AC network and the results could be compared with the DC circuit. Another direction is to expand the model by incorporating fuel-based DG or batteries and analyze the tradeoff between the cost and the resilience performance in various operating scenarios.

APPENDIX
Loss of load probability (LOLP) is defined as the probability that the power supply drops below the load. We rewrite the chance constraint (15) as follows Pr{P t (x) < D t } ≤ α 1 , ∀t. (A1) where, α 1 is the acceptable LOLP limit. Define Y = P t (x) − D t , then we have Pr{Y < 0} ≤ α 1 , ∀t.
(A2) and the mean and variance of Y is According to central limit theorem, P t (x) and D t follows the normal distribution, so does Y . Therefore, the deterministic counterpart of constraint (A1) or (A2) is obtained as where Z 1−α 1 is the z-score at (1-α 1 )×100% confidence for a standard normal distribution. By substituting Equations (A3) and (A4) into (A5), we have Similarly, under the normality condition the deterministic constraints (25) and (26) for voltage variation and thermal stress control are derived based on the corresponding chance constraints (16) and (17), respectively.