A Novel Decomposition Solution Approach for the Restoration Problem in Distribution Networks
Abstract β The distribution network restoration problem is by nature a mixed integer and non-linear optimization problem due to the switching decisions and Optimal Power Flow (OPF) constraints, respectively. The link between these two parts involves logical implications modelled through big-M coefficients. The presence of these coefficients makes the relaxation of the mixed-integer problem using branch-and-bound method very poor in terms of computation burden. Moreover, this link inhibits the use of classical Benders algorithm in decomposing the problem because the resulting cuts will still depend on the big-M coefficients. In this paper, a novel decomposition approach is proposed for the restoration problem named
Modified Combinatorial Benders (MCB) . In this regard, the reconfiguration problem and the OPF problem are decomposed into master and sub problems, which are solved through successive iterations. In the case of a large outage area, the numerical results show that the MCB provides, within a short time (after a few iterations), a restoration solution with a quality that is close to the proven optimality when it can be exhibited.
Index Terms β Convex Optimization Problem, Decomposition, Distribution Network, Load Pickup, Line Switches, Reconfiguration, Restoration Service.
NOMENCLATURE A. Parameters π€ ππ , π€ π π€ , π€ ππ Weighting factors of the objective function terms (p.u.) π· π Importance factor of the load at bus π (p.u.) π ππ (π π ) The operation time of line switch ππ /load breaker π (hour). π π,π‘π· (π π,π‘π· ) Active (/Reactive) power demand at bus π, at time π‘ (p.u.). π ππ (π₯ ππ ) Resistance (Reactance) of line ππ (p.u.). π£ πππ₯ (π£ πππ ) Maximum (Minimum) limits of voltage magnitude (p.u) π πππππ₯ Maximum current flow rating of line ππ (p.u) π π,πππ₯πππ Active power capacity of resource at node π (p.u.) π π,πππ₯πππ Apparent power capacity of resource at node π (p.u.) π A large multiplier B. Variables π ππ Binary decision variable indicating if the line ππ equipped with a switch is energized or not (1/0) πΌ π,π‘ Binary decision variable indicating if at time π‘ the load at node π is supplied or rejected (1/0) π½ ππ Continuous variable indicating if the line ππ is oriented from node π to node π or not. π ππ Auxiliary flow that is travelling in line ππ from node π to node π . π π Auxiliary variable indicating if the node π is energized or not (1/0). πΉ ππ,π‘ Square of current flow magnitude in line ππ, at time π‘ (p.u) π ππ,π‘ (π ππ,π‘ ) Active (Reactive) power flow in line ππ , starting from node π , at time π‘ (p.u). π π,π‘πππ (π π,π‘πππ ) Active (Reactive) power injection from the substation or DGs at node π , at time π‘ (p.u.). π π,π‘ Square of voltage magnitude at bus π, at time π‘ (p.u). C. Indices π, π
Index of nodes ππ Index of branches π‘ Index of time π Index of iteration in the decomposition approach D. Sets π Set of nodes π β Set of nodes in the off-outage area π Set of lines (plus tie-lines) π β Set of lines (plus tie-lines) in the off-outage area π π =π ππ£ππ βͺπ πππ‘π βͺπ π πππ Set of lines (plus tie-lines) in the off-outage area equipped with switches π ππ£ππ Set of lines hosting available tie-switches π πππ‘π Set of lines hosting internal tie-switches π π πππ Set of lines hosting sectionalizing switches πΊ π·πΊ Set of nodes hosting DGs I. I NTRODUCTION
After a failure in a radial distribution network, once the fault is isolated, the area downstream to the fault place remains unsupplied. This area is called off-outage area . The aim of the restoration operation is to restore the maximum energy of loads within this off-outage area while minimizing the total switching operation time [1]. In order to achieve this goal in passive distribution networks, the only possible action is to transfer the unsupplied loads to the healthy neighboring feeders (Fig. 1.). This reconfiguration is deployed through changing the status of normally-closed (sectionalizing) and normally-open (tie) switches. The tie-switches that are between the faulted feeder and a healthy feeder are called available tie-switches (T2, T3, and T4 in Fig. 1.). The tie-switches with both ending nodes inside the off-outage area are referred to as internal tie-switches (T5 in Fig. 1.) . The resulting new configuration of the network remains for a so-called restorative period that starts from the fault isolation instant until the time when the faulted element is repaired. After the restorative period, the original configuration of the network will be restored. The highly increasing penetration of Distributed Generators (DGs) in Active Distribution Networks (ADNs) introduces, among others, new restoration actions besides the switching operations. Among these additional restoration actions are DG power set point modifications. The incorporation of these actions could lead to a more efficient restoration solution.
A Novel Decomposition Solution Approach for the Restoration Problem in Distribution Networks
Hossein Sekhavatmanesh,
Student Member, IEEE , Rachid Cherkaoui,
Senior
Member, IEEE T4 F ee d er - a F e e d e r - c
146 7 8 9 5
Feeder-b
T1T2 T3
15 16 T5 Circuit breakerTie-switchBus or load center
Fault
Sub2Sub1 DG Cluster 1Cluster 2
Fig. 1. A simple distribution network under post fault conditions.
However, this leads to increases the complexity of the restoration problem. Therefore, the challenge is to find an automatic and efficient solution for the restoration problem as an operator decision support in distribution networks while considering their new active status. From mathematical formulation point of view, the restoration problem is an NP-hard combinatorial optimization problem. It is in the form of a mixed-integer and non-linear problem, respectively, due to the switching decisions and Optimal Power Flow (OPF) constraints. The OPF is a known challenging optimization problem. It has been the main building block for the formulation of many control, operation, and planning problems in power systems. In this respect, an AC-OPF should be integrated into the restoration problem in order to check the feasibility of the obtained solution concerning the technical constraints such as voltage and current limits. In order to deal with the non-polynomial hardness of such an OPF-based optimization problem, some researchers applied meta-heuristic methods as the solution approach. Among these methods are Genetic Algorithm (GA) [2], Particle Swarm Optimization (PSO) [3], Ant Colony [4], and Tabu Search [5]. The main problem with these methods is that they are in general computationally exhaustive and could fail to find an existing feasible and good-enough solution in a reasonable time complying with online restoration requirements. In order to address this weakness, some researchers proposed to use heuristic approaches. These approaches use graph search methods to find a suitable topology of the network with regard to the restoration criteria [6], [7]. It means that the heuristic restoration algorithms will be applicable only for specific network topologies and cannot be generalized for any network topology. Moreover, as it is shown in [8], the restoration solution provided by these heuristic algorithms could be very far from the global optimal solution. The mathematical analytical programming started to be an option for solving the restoration problem shortly after that some convex relaxation methods were proposed for the OPF problem. The authors of [9]-[10] used semidefinite programming (SDP) relaxation in order to solve it for radial unbalanced grids. Regarding the grid unbalancing, a distributed optimization technique is developed in [11] based on the Alternating Direction Method of Multipliers (ADMM). In this paper, we assume that the distribution network is operated in a radial and balanced fashion. In this respect, relaxation methods are proposed in the literature for the OPF problem in the form of either Mixed-Integer Quadratic Constraint Programming (MIQCP) [12], [13] or Mixed-Integer Second-Order Cone Programming (MISOCP) model [14]β[16]. The discussion given in Appendix I gives an insight about possible extension of the proposed restoration approach to unbalanced grids. The convexity of the power flow formulation allows finding a solution with proven optimality up to the desired accuracy. However, the computation burden of the resulting optimization problem could be intractable depending on the dimension of the distribution network. This drawback inhibits the use of mathematical programming methods to tackle the restoration problem as a multi-period and/or multi-scenario optimization problem. In order to address this problem, different studies proposed decomposition methods. Some of them account for the uncertainty of load demands and DG injections using stochastic or robust optimization [17]β[21]. They proposed decomposition algorithms where the innermost problem is assigned to find the worst realization of uncertain parameters given a fixed radial configuration. In the uppermost level, the deterministic restoration problem is solved while fixing the uncertain parameters to their worst-case realizations found in the inner level problem. For solving this decomposed restoration problem, different solution approaches have been applied. Among the most important ones are stochastic rolling-horizon optimization method [18], Information Gap Decision Theory (IGDT) [19], and column-and-constraint generation algorithm [20], [21]. The restoration strategy should consider the time-varying loads in order to provide a unique configuration valid throughout the whole restorative period. Actually, none of the above-mentioned papers accounts for the time-varying loads. If those decomposition algorithms are used for solving a such multi-period optimization problem, the deterministic restoration problem in the uppermost level will be computationally intractable. In this regard, papers [22], [23] propose to partition the time window of the optimization problem in clusters with similar load levels. Then, the reconfiguration problem is solved sequentially for each cluster of time instants. The weakness of this strategy is that the solutions at a given sequence are not influenced by the future sequences during the rest of the optimization process. Therefore, this approach cannot be applied to solve multi-period restoration problems (such as the one studied in this paper), where the feasible and optimal solution at one time step depends on the solution of the problem in the previous and next steps. In order to relax the computation burden of the multi-period restoration problem, the authors in [17], [24], [25] propose to use a linear load flow formulation instead of the original and non-linear formulation. Since in the restoration strategy, the reconfigured network could be operated very close to its capacity envelop, applying the linear load flow model may result in an infeasible solution regarding the network safety constraints (e.g. voltage and current limits). Another approach for relaxing the computation burden of the restoration problem is provided in [26]. According to this methodology, the restoration problem is solved in two stages. In the first stage, the post-restoration topology is determined using a heuristic approach. The set of loads to be restored and the outputs of sources are determined in the second stage while fixing the network topology to the one obtained in the first stage. The optimization problems in the first and second stages, which are referred as reconfiguration problem and load pickup problem, respectively, are mutually interdependent. The decoupling of these two problems as proposed in [26] could lead either to no feasible solution or to a solution very far from the optimal one. In order to address all the afore-mentioned weaknesses, we propose a modification to the combinatorial Benders method so that it can be used for the multi-period restoration problem while considering case studies of realistic size. Combinatorial Benders method was firstly proposed by Hooker for solving optimization problems which include conditional constraints [27, p.]. This method has been used in many different applications such as in circuit verification problems [28], map labeling problem [29] and asymmetric travelling salesman problem [30]. In this paper, a two-stage analytical formulation is proposed for the restoration problem. It is solved thanks to a so called Modified Combinatorial Benders algorithm. Compared with the state-of-the-art, the major contributions of this paper are the following: 1-
A novel decomposition approach is proposed for the restoration problem while considering inter-temporal constraints (e.g. varying loads) and control actions (e.g. DG power set points). The original problem is in the form of a multi-period, mixed integer, and non-linear optimization problem. Thanks to the proposed decomposition approach, the restoration problem is made tractable for analytical solvers in case of a grid of realistic size in a multi-period optimization problem. 2-
The standard combinatorial Benders method is augmented with new cuts identifying binary variable combinations that are either infeasible or non-optimal. In this regard, the proposed cuts distil the search space of the optimization problem at a given iteration into a smaller subset that includes the global optimal solution. Therefore, compared with the standard combinatorial Benders, the proposed MCB approach converges in less number of iterations. The quality of this final solution is close to the global optimal solution. 3-
A convex AC-power flow formulation is integrated into the decomposed formulation proposed for the restoration problem in order to accurately model the electrical operational constraints (e.g. voltage and current limits). The remainder of this paper is organized as follows. The integrated mathematical model of the restoration problem is presented in section II. The decomposition approach is presented in section III. In this section, first, the details of the Combinatorial Bender decomposition method are explained. Then, the proposed modified one is provided. Section IV illustrates different case studies verifying the main advantages of the proposed decomposition algorithm with respect to the integrated analytical optimization method. Finally, section V concludes the paper highlighting the main contributions. II. P ROBLEM F ORMULATION
In this section, an integrated model of the time-dependent restoration formulation is presented in the form of a MISOCP problem. This problem encompasses three groups of decision variables namely, I) the binary variables π ππ represent the energization status of line switch ππ , II) the binary variables πΌ π,π‘ account for the status of the load at node π at time π‘ , III) and the continuous variables π π,π‘π·πΊ /π π,π‘π·πΊ are associated to the active/reactive power set points of DG at node π at time π‘ . The targeted restoration strategy is a multi-period optimization problem in the sense that the decision for the load pickup and DG power set points varies with time. The line switching variables ( π ππ ) do not vary with time, since it is aimed to provide a single new network configuration for the whole restorative period. Minimize: πΉ πππ = π ππ . πΉ ππ + π π π€ . πΉ π π€ + π ππ . πΉ ππ (1.a) where, πΉ ππ = β β π· π . (1 β πΌ π,π‘ ). π π,π‘π· πβππ‘ πΉ π π€ = β π ππ . π ππ(π,π)βπ π‘πππ + β (1 β π ππ ). π ππ(π,π)βπ π πππ πΉ ππ = β β π ππ . πΉ ππ,π‘(π,π)βππ‘ Subject to: { 0 β€ π½ ππ β€ 1 β(π, π) β π S π½ ππ + π½ ππ = π ππ , β(π, π) β π π π½ ππ = π ππ , π½ ππ = 0, β(π, π) β π ππ£ππ (1.b) {π π = β π½ πππ:(π,π)βπ β β€ 1, βπ β π β π π = 1, βπ β π\π β (1.c) { 0 β€ Ξ¨ ππ β€ π. π½ ππ β(π, π) β π S ππ β€ π. π½ ππ β(π, π) β π S β (Ξ¨ π β π ) βπ β :(π β ,π)βπ π = β (Ξ¨ ππ β ) βπ β :(π,π β )βπ π + π π βπ β π β β Ξ¨ ππβ(π,π)βπ ππ£ππ = β π πβπβπ β (1.d) {πΌ π,π‘ β€ π π , βπ β π β , βπ‘ πΌ π,π‘ = 1, βπ β π\π β , βπ‘ (1.e) π,π‘β1 β€ πΌ π,π‘ β€ 1 βπ β π\π β , βπ‘ (1.f) βπ. (1 β π ππ ) β€ π π,π‘ β π π,π‘ β 2(π ππ . π ππ,π‘ + π₯ ππ . π ππ,π‘ ) β€ π. (1 β π ππ ) βππ β π, βπ‘ (1.g) { π ππ,π‘ = ( β π ππ β ,π‘π β β π(π β ,π)βπ )+π ππ . πΉ ππ,π‘ + πΌ π,π‘ . π π,π‘π· β π π,π‘πππ π ππ,π‘ = ( β π ππ β ,π‘π β β π(π β ,π)βπ )+π₯ ππ . πΉ ππ,π‘ + πΌ π,π‘ . π π,π‘π· β π π,π‘πππ (1.h) βππ β π, βπ‘ β 2π ππ,π‘ ππ,π‘ πΉ ππ,π‘ β π π,π‘ β β€ πΉ ππ,π‘ + π π,π‘ β(π, π) β π, βπ‘ (1.i) { 0 β€ π π,π‘πππ β€ π π,πππ₯πππ βπ π,π‘πππ π π,π‘πππ β β€ π π,πππ₯πππ βπ β Ξ© π·πΊ β© π, βπ‘ (1.j) { 0 β€ πΉ ππ,π‘ β€ π ππ . π πππππ₯2 βπ. π ππ β€ π ππ,π‘ β€ π. π ππ βπ. π ππ β€ π ππ,π‘ β€ π. π ππ βππ β π, βπ‘ (1.k) π£ πππ2 . π π β€ π π,π‘ β€ π£ πππ₯2 . π π βπ β π, βπ‘ (1.l) The objective function (1.a) tends to minimize the weighted total costs associated with the reliability ( πΉ ππ ), switching ( πΉ π π€ ) , and operational ( πΉ ππ ) objectives, in decreasing order of priority. This hierarchical priority is enabled using π β ππππ π‘πππππ‘ method [31]. The reliability cost is the total energy not supplied of the loads while accounting for their importance factors. The switching cost is formulated as the summation of two sub terms associated with the total operation time of tie-switches and sectionalizing switches, respectively. Finally, the least priority objective term is the operational term formulated as the total active power loss. As the formulation of the restoration problem is not the focus of this paper, only a very brief description of constraints (1) will be provided in the following. The readers are recommended to go through [8] for detailed explanation. Constraints (1.b)-(1.d) model the reconfiguration problem ensuring the radial topology of the network [32]. Constraints (1.d) ensure that all the nodes are connected to the substation via flows of a virtual commodity represented by auxiliary and continuous variables Ξ¨ ππ and Ξ¨ ππ . A brief explanation of these constraints is provided in Appendix IV. Using these constraints, we avoid to create an isolated loop in the off-outage area with the loads that are supplied in an islanded way using an existing DG in that loop. Constraints (1.e)-(1.f) formulate the load pickup problem. For an energized node π in the off-outage area ( π π = 1 ), a decision is made in (1.d) with binary variable πΌ π,π‘ , indicating if its load is restored at time t or rejected (1/0). As formulated in (1.f), it is assumed that once a given load is restored at a given time, no further interruption is permitted during the subsequent time slots of the restorative period. The relaxed formulation of AC-OPF is presented in (1.g)-(1.j). The aim of this part is to dispatch the active/reactive power set points of DGs, while respecting all the security constraints in the reconfigured network. The reconfiguration problem is linked to the AC-OPF problems in (1.k) and (1.l) using variables π ππ and π π , respectively. These links inhibit the use of classical Benders algorithm in decomposing the problem because the resulting cuts will still depend on the big-M coefficients used in (1.k) and (1.l). In the subsequent discussion, it is aimed to present a tractable approach for solving the restoration problem in a multi-period optimization environment. In this respect, the following compact form of the restoration model will be used to represent the above extensive formulation. For example, if line ij is energized ( π ππ = 1 ), then the current flow in this line must be less than its ampacity limit ( πΉ ππ,π‘ β€ π πππππ₯2 ). π β min π₯,π¦,π§,π’ β π π (1 β π₯ π ) πβπ β + β π ππ π¦ ππππβπ β + β π ππ π§ ππππβπ (2.a) Subject to: β (π₯, π¦) β β (2.b) β (π₯, π¦, π§, π’) β β (2.c) where: β (π₯, π¦) β { π¦ β Ξπ΄π₯ β₯ ππ΅π¦ + πΆπ₯ β₯ π (3.a) (3.b) (3.c) β (π₯, π¦, π§, π’) β { ππ π π (π¦) = 1 π‘βππ π· π π§ β₯ π π , βπ = 1,2, . . . , π πππ₯ πΈπ§ = π½π₯ β πΉπ’ βπΊ π π§β β€ π ππ π§, βπ β π βπ» π π’β β€ β π , βπ β Ξ© π·πΊ (4.a) (4.b) (4.c) (4.d) where, π¦ and π₯ are vectors of binary variables indicating, respectively, line switching variables ( π ππ ) and load pickup variables ( πΌ π,π‘ ) . Continuous variables are represented by vectors π’ and π§ , which refer, respectively, to the DG power set point variables ( π π,π‘π·πΊ /π π,π‘π·πΊ ) , and the rest of state variables related to the optimal power flow constraints at each time t (such as π π,π‘ , πΉ ππ,π‘ , β¦ ). The three terms of (2.a) represent, respectively, the reliability, switching, and operational objective terms formulated in (1.a). β is expressed in (3) as the set of constraints only on the binary variables ( π₯ and y). In (3.a), Ξ is the set of radial network configurations described by (1.b)-(1.d). Constraint (3.b) accounts for (1.f) as the load pickup formulation. Constraint (3.c) represents (1.e) as the link between reconfiguration and load pickup problems. As given in (4), β represents the set of AC-OPF constraints which are linked to the binary variables (x and y). The link between the reconfiguration and AC-OPF problem is given in (4.a) in the form of conditional constraints, which are linearized using big-M formulation as formulated in (1.k) and (1.l). They mean that if a certain condition on π¦ variables holds ( π π (π¦) =1 ), a constraint on π§ variables is added to the optimization problem . Equation (4.b) accounts for the voltage equation and the power balance equation given in (1.g) and (1.h), respectively. Equations (4.c) and (4.d) represent second-order constraints associated with the current flow equation (1.i), and DG apparent power capacity limit (1.j). In (4.d), β π refers to the apparent power capacity of the DG at node π . The other notations used in (2),(3) and (4) are matrices (the ones in capital letter) or vectors (the ones in small letter) of parameters. III. S OLUTION S TRATEGY
In this section, a novel decomposition approach is proposed for the restoration problem named Modified Combinatorial Benders. In this regard, the reconfiguration problem and the OPF problem are decomposed into master and sub problems, which are solved through successive iterations. A. Master Problem
In the outer level of the proposed decomposition strategy, the master problem is solved. The master problem is the same as the original problem π while removing the AC-OPF constraints ( β ). However, the operational constraints (e.g. voltage and current limits) are not completely disregarded as they are represented by DistFlow constraints (β Μ Μ Μ ). β³ (π) β min π₯,π¦,π§,π’ β π π (1 β π₯ π ) πβπ β + β π ππ π¦ ππππβπ β (5.a) Subject to: β (π₯, π¦) β β (5.b) β Μ Μ Μ (π₯, π¦, π§, π’) β β (5.c) { β π π (1 β π₯ π ) πβπ(π,π¦ (π) ) β₯ β Ξ½ (π¦ (π) ) βΆ π¦ β π(π, π¦ (π) ), β Ξ½ (π¦ (π) ) β β π¦ β π(π, π¦ (π) ) βΆ β Ξ½ (π¦ (π) ) = β βπ = 1,2, β¦ , π πππ₯ (π) , π = 1, β¦ , π β 1 (5.d) where: β Μ Μ Μ (π₯, π¦, π§, π’) β { ππ π π (π¦) = 1 π‘βππ π· π π§ β₯ π π , βπ = 1,2, . . . , π πππ₯ πΈπ§ = π½π₯ β πΉπ’ πΊ π Μ π π§ β€ π Μ π , βπ β Ξ© ππ’π π»Μ ππ π’ β€ βΜ π , βπ β Ξ© π·πΊ (6.a) (6.b) (6.c) (6.d) β Μ Μ Μ is formulated in (6) as the set of linear DistFlow constraints linked to the binary variables. Constraints (6.c) and (6.d) are the linearized formulation of line ampacity and DG apparent power capacity limits, where π π Μ and β π Μ induce relaxed line ampacity limit and DG apparent power limit, respectively. This linearization technique is according to the technique presented in [20]. The detailed formulation of the DistFlow constraints is given in [33]. Two sets of constraints (5.d) are denominated as optimality cuts and feasibility cuts, respectively. These constraints represent the modified version of the Combinatorial Benders cuts introduced in [34]. The main idea behind this modification is to have recourse functions providing information not only on the feasibility of each solution but also on its optimality. According to the standard Combinatorial Benders method, all the binary variables must be set as the complicating variables. They are determined by the master problem and the sub problem just to evaluate the feasibility of the solution. In case that the solution is infeasible, a cut will be added to the master problem for the next iteration to remove the corresponding set of infeasible binary variables from the solution space. In the proposed modified version, only a subset of binary variables is fixed in the sub problem and defined as the complicating variables. The other binary variables that are not fixed in the sub problem are called floating variables. In the formulation provided in (5) π¦ and x represent, respectively, complicating and floating variables. In the case where the solution of the master problem at iteration k ( π¦ (π) ) leads to no feasible solution in the sub problem, the feasibility cuts are augmented by the second constraint of (5.d). If the sub problem at iteration k is feasible, its optimal solution β Ξ½ (π¦ (π) ) is used to augment the optimality cuts according to the first constraint of (5.d). The constraints given in (5.d) are non-convex and need to be linearized. This linearization together with the formulation of π(π, π¦ (π) ) is provided in section III.C. πΏπ΅ (π) = β³ (π) (7) As given in (7), the master problem β³ (π) formulated in (5.c) provides a lower bound for the original optimization problem π expressed in (2). Actually, unlike in the case of AC-OPF formulation in the sub problem, we choose relaxed limits for the electrical operation constraints (e.g. voltage limits) of the DistFlow formulation in the master problem. Therefore, the feasible region of π (π) under DistFlow constraints is relaxed in comparison to the feasible region of π· under AC-OPF constraints. B. Sub Problem
Once the optimal configuration is found in the master problem, the next step is to find the optimal load pickup solution, if any, for the obtained configuration subject to the AC-OPF constraints. When we fix the network configuration, the topology of the off-outage area will be partitioned accordingly into several clusters. A cluster is defined as a collection of nodes and lines in the off-outage area that are supplied by only one available feeder.
π(π, π¦ (π) ) determines the index of nodes that are in cluster π . π(π, π¦ (π) ) denotes a set of configurations which provide the same optimal load pickup solutions for the nodes in cluster π . The formulation of π , optimality cuts and feasibility cuts are provided in section III.C. Consider the simple network of Fig. 1 as an example. The switching operations shown in this figure are assumed to represent the optimal solution found in the master problem. Under the resulting configuration, the off-outage area is partitioned in two clusters. First cluster includes nodes 14, 15 and 16 that are restored from feeder-c through tie-switch T4. The second cluster includes nodes 11, 12 and 13 that is supplied from feeder-b through tie-switches T2 and T3. As shown in Fig. 1., there is no path between two nodes belonging to the two different clusters except through the slack bus. We assume that the slack buses are effectively fixing the voltage set point at the top of each feeder (bus βSub2β in Fig. 1.). Under this assumption, it can be said that the change of loading in one cluster (change of π₯ variables) does not change any state variable outside that cluster. Therefore, we solve a separate sub problem for each individual cluster. The aim is to break the computation burden of the inner level problem into several problems, which can be handled using different cores in parallel. The following MISOCP formulation models the sub problem for cluster π at iteration q , given network configuration π¦ (π) . β Ξ½ (y (q) ) β min π₯,π§,π’ β π π (1 β π₯ π ) πβ π ( π,π¦ (π) ) + β π ππ π§ ππππβ π(π,π¦ (π) ) (8.a) Subject to: β (π₯, π¦ (π) ) β β (8.b) β (π₯, π¦ (π) , π§, π’) β β (8.c) where, π denotes the index of lines within cluster π . It should be noted that the sub problem β π(q) incorporates only those variables that are related to the lines and nodes in cluster π . The objective function (8.a) is the minimization of the total energy not supplied in cluster π as formulated in (1.a). Since the complicating variables π¦ are fixed, the big-M coefficients used in (4.a) do not appear in the sub problem which is relaxing again the computation burden of the inner level problem with respect to the original optimization problem π given in (2 ). ππ΅ (π) = β β π (π¦ (π) ) π πππ₯ π=1 (9) According to (9), the summation of optimal objective values β Ξ½ (y (q) ) , associated with all the clusters, induce an upper bound ππ΅ (π) for the reliability optimal solution in the original optimization problem π (2). The optimal solution of the sub problem β Ξ½ (y (q) ) is also used to augment the feasibility cuts of the master problem as formulated in the first constraint of (5.d). In case that there is no feasible solution for the sub problem, a feasibility cut is generated and added to the master problem as formulated in the second constraint of (5.d). C. Generating the optimality and feasibility cuts
In a given iteration k, if the sub problem β Ξ½ (π¦ (π) ) has a solution, say π₯ (π) , then clearly, π₯ = π₯ (π) is the optimal solution of P if y= π¦ (π) . We look for π as a set of y-solutions, such that if we solve the sub problem while fixing y variable to any point π¦β in this set, no better solution than π₯ (π) can be found. It means that for any yβ² β π , β Ξ½ (π¦β²) β₯ β Ξ½ (π¦ (π) ) . This constraint is formulated in the first expression of (5.d) and denominated as an optimality cut. If the sub problem β Ξ½ (π¦ (π) ) at iteration k is infeasible, then π is defined as a set of π¦ -solutions, such that if π¦ variable is fixed to any other point π¦β in this set, the sub problem will be still infeasible. In other words, in order to break the infeasibility, variable π¦ should take values outside the set of π . This constraint is formulated in the second expression of (5.d) and referred as the feasibility cut. Note that π(π, π¦ (π) ) is associated with a given master problem solution π¦ (π) and also with a cluster π . Also note that the solution of the master problem, say π¦ (π) , represents the network configuration and the solution of sub problem, say π₯ (π) , is the value of load pickup variables in cluster π . According to the definition of the optimality and feasibility cuts, in order to derive π , we should find network configurations π¦β that lead to no better reliability solutions in cluster π , with respect to β Ξ½ (π¦ (π) ) . The optimal solution of load pickup variables within cluster π will not improve under configuration π¦β with respect to the optimal values under configuration π¦ if the following conditions hold: a) All the nodes in cluster π that were connected to each other under configuration π¦ (identified by π (π, π¦ π) ) ), are still connected to each other under configuration π¦β . b) The injection nodes that were supplying the nodes in cluster π under configuration π¦ are the same as those under configuration π¦β² , Consider the test system of Fig. 1, as an example. As mentioned earlier, nodes 14, 15 and 16 are in the first cluster that is supplied by feeder-c through tie-switch T4 and by DG at node 5 through tie-switch T5. Tie-switches T4 and T5 are named as source lines. Source lines of cluster π are defined as the lines at the border of cluster π that are injecting active or-/and reactive power to the cluster. Considering the example shown in Fig. 1, assume that all the nodes in the first cluster should be restored except node 16, according to the solution of the sub problem β Ξ½ (π¦ (π) ) . Now, by changing the configuration, it is obvious that the load at node 16 still cannot be restored if a) nodes 14, 15 and 16 are still connected to each other and if b) these nodes are supplied through the same source lines (tie-switches T4 and T5). According to two conditions mentioned above, the set π (π, π¦ (π) ) is expressed as in (10). π( π, π¦ (π) ) = {π¦|βπ β² β€ π πππ₯ ( π¦ ) βΆ π (π, π¦ (π) ) β π (πβ², π¦), Ξ₯(π β² , π¦) β Ξ₯(π, π¦ (π) )} (10) where, Ξ₯ represent the index of source lines that are injecting power to cluster π . In order to preserve the linearity of the optimality and feasibility cuts in terms of π¦ variables, the two constraints expressed in (10) are reformulated in (11) and (12). β π¦ ππβπ(π,π¦ π) ) = | π (π, π¦ (π) )| β 1 (11) π¦ π β€ π¦ π(π) βΆ βπ β π(π, π¦ (π) ) (12) The connectivity condition mentioned in condition a) is formulated in (11). This constraint enforces that the number of closed lines in a given cluster π is equal to the total number of nodes in cluster π minus one. This is the tree condition for cluster π . The tree condition ensures the network connectivity if it is radial [8]. This radiality condition is ensured for a given cluster π using (5.b). Constraint (12) formulates condition b) that is mentioned above. This constraint ensures that the resource line π of cluster π that was open under configuration π¦ (π) will stay open under any configuration π¦ β π (π, π¦ (π) ) . According to the derived formulations for the set of π (π, π¦ (π) ) , the feasibility cut that was given in the second constraint of (5.d) is re-formulated in the following. β π¦ ππβπ(π,π¦ π) ) β€ |π(π, π¦ (π) )| β 2 (13) β π¦ ππβπ(π,π¦ π) ) β₯ 1 + β π¦ πππβπ(π,π¦ π) ) (14) where, at least one of the conditions (13) or (14) must hold. This either-or constraint cannot be integrated into a convex model. Since in a convex optimization problem, all the constraints must hold. Therefore, this constraint should be further re-formulated according to the strategy given in Appendix II. Expressions (13) and (14) are the complements of (11) and (12), respectively. It should be mentioned that the complement of (12) means that at least one line π β π(π, π¦ (π) ) exists such that π¦ π β₯ π¦ π(π) + 1 . This constraint can be expresses as in (14). Regarding the optimality cut, the first expression given in (5.d) can be translated into the following: If (11) and (12) are satisfied, then β π π (1 β π₯ π ) πβπ(π,π¦ (π) ) β₯ β π (π¦ (π) ) must be satisfied. Therefore, the optimality cut is in the form of a conditional constraint. In order to be integrated in a convex optimization model, it is re-formulation as given in Appendix III. D. Modified Combinatorial Benders Algorithm
The proposed decomposition approach for solving the distribution network restoration problem is described as follows: 1-
Initialize iteration number ( π β 1 ), lower bound (
πΏπ΅ β 0 ), upper bound (
ππ΅ β +β ), and set the convergence tolerance ( π > 0 ). 2-
Solve the master problem to get the optimal network configuration π¦ (π) . Update the lower bound ( πΏπ΅ βmax (πΏπ΅, πΏπ΅ (π) ) ), where πΏπ΅ (π) is given in (7 ). 3- Solve the sub problem for the obtained configuration π¦ (π) and for each cluster π . a. If the optimization problem is feasible, find the optimal load pickup variables π₯ (π) and augment the optimality cuts according to the first constraint of (5 .d). Update the upper bound ( ππ΅ β min (ππ΅, ππ΅ (π) ) ), where ππ΅ (π) is given in (9 ). b. If the optimization problem leads to no feasible solution, augment the feasibility cuts according to the second constraint of (5 .d). 4- Check for convergence : a. If ππ΅ β πΏπ΅ β€ π πππ‘ or if the computation time is larger than π π‘πππ , then terminate the algorithm and propose the best UB solution found so far as the solution of the problem. b. Else, update the iteration number ( π β π + 1 ), and return to step 2. While the iterations of the proposed algorithm are evolving, the original solution space is gradually reduced by removing more combinations of binary variables. This is realized in a conservative way using the proposed optimality and feasibility cuts. Therefore, using the MCB approach, we might not be able to converge to the global optimal solution. However, as it will be illustrated in section IV.A, when the MCB algorithm converges, the best solution visited so far is close to the global optimal solution. In order to end up with a tractable solution methodology in case of grids with realistic sizes, two stopping criteria are defined in the above-mentioned algorithm. According to this algorithm, we continue the running of iterations until the difference between the lower bound solution (LB) and the upper bound solution (UB) is lower than a threshold ( π πππ‘ ). In addition, we impose an additional threshold on the computation time ( π π‘πππ ). In this regard, if the computation time is more than a threshold value, then the algorithm is stopped. The values of these thresholds ( π πππ‘ and π π‘πππ ) are determined based on the experience of DSO. IV. N UMERICAL A NALYSIS AND D ISCUSSION
In order to illustrate different features of the proposed solution algorithm for the restoration problem, two medium voltage networks are used shown in Fig. 2, and Fig. 3. In this paper, we study different scenarios. Scenarios I and II are applied on the test network of Fig. 2, whereas for scenarios III,
F1F3
WTDGPV DG OLTC2OLTC1 DG T1 B1 F ee d er A B2 F e e d e r B
15 20 2218 24 B3 Feeder C
25 29 B4 F ee d er D
31 34 B5 F ee d er E B6 F ee d er F
502 60247 B7 F ee d er G B8
64 6570 B9 F ee d er H F ee d er I B10 F e e d e r J B11 T2T3
To 13
T4 T5T6 T8T9T10 T12 F e e d e r K T11T7 DG B12 B13 B15B14 F2 Sectionalizing switchCircuit BreakerTie-switchBus or load centerCritical load F6 Fig. 2. The test network for test scenarios I and II [4].
Feeder1 Feeder2Feeder3 Feeder4
T5T7T2 T1T4 T8 T3
B1 B2B3 B414
T6 F4
Sectionalizing switchCircuit BreakerTie-switchBus or load centerCritical loadDispatchable DG z zz F5 Fig. 3. The test network for test scenario III. the test network of Fig. 3 is used. For both test networks, the base power and energy values are assumed to be 1MW, and 1MWh, respectively. The minimum and maximum voltage magnitude limits are set, respectively, to 0.917 and 1.05 p.u [35]. The hourly profiles for different types of load patterns that are used in the both test networks are given in [36]. Two types of DGs are considered in this paper. Namely dispatchable and non-dispatchable. The dispatchable DGs, such as the diesel generators, are controlled to deliver the active and reactive power references that are set by DNO ahead of their operation. We consider also non-dispatchable DGs such as PV and wind generators, which are modeled as voltage-independent active power injection units with zero reactive power components. The forecast power injections of PV- and Wind-based DGs are derived from the real data reported in [37] and [38], respectively. In both test network, it is assumed that each node is equipped with a load breaker. All the line switches are
Table I. Comparison of restoration re sults obtained using IAO and MCB methods in scenarios I and II. S ce n a r i o S o l u ti on M e t hod R ec on f i gu r a ti on A c ti on s L o a d P i c kup A c ti on s πΉ π π ( p . u . ) πΉ π π€ ( m i n ) C o m pu t a ti on ti m e ( s ec ) I (f a u lt s F a nd F ) I AO I. Open switch 38-39 and load breakers {33,34,37,41,42}, and close T11 II. Close T3 - 88 92.5 2.12 M CB I. Open switch 35-36 and load breakers {33,34,37,38,41,42}, and close T11 II. Close T3 - 100.1 93 9.58 II (f a u lt F ) I AO I. Open switches {1-2 ,33-34,31-41} and load breakers {2,3,4,5,6,12,13,14,35,36,37,38,39, 40,41}, and close {T3,T11} II. Close{T1,T4} III. Close load breaker 37 at t=12 IV. Close load breaker 13 at t=19 5.97e3 219 385.3 M CB I. Open switches {6-7 ,11-12,31-41} and load breakers {1,3,12,14,15,17,19,20,22,25,26,27, 28,29, 30,31,32,34,35,36,38, 43,44,46}, and close {T3,T7,T10} II.Close{T1,T2,T5,T8} III. Close load breakers {3,22,25,29} at t=19 IV. Close load breaker 13 at t=20 4.13e3 313.8 120 assumed manually-controlled, whereas the load breakers are all assumed remotely-controlled. The time needed for the operation of each manually controlled and remotely-controlled switch are assumed 30 and 0.5 minutes, respectively. It is assumed that the critical loads that are shown with β β β. in Fig. 2 and Fig. 3, have the priority factors equal to 100 while the priority factors of other loads are equal to 1. In order to show the functionality of the proposed solution approach, the restoration problem in case of each case study is solved using two approaches: I) the Integrated Analytical Optimization (IAO) method, and II) the Modified Combinatorial Benders (MCB) decomposition method proposed in section III. According to the IAO approach, the integrated optimization problem (1) is solved in one shot using an analytical solver. For this aim, the Branch-and-Bound method is used to relax the integrality constraints of the original optimization problem in an iterative way. In this regard, the best integer (valid) solution that is found at any step in the algorithm is called incumbent solution. The objective value of this incumbent is an upper bound for the optimal solution of the original minimization problem. At any step through the Branch-and-Bound search algorithm, there is also a lower bound, called the best-bound solution. This bound is obtained by taking the minimum of the optimal objective values of all the solutions obtained so far including the infeasible ones regarding the integrality constraints. The difference between the current upper Fig. 4. The progress of the obtained solution using IAO and the MCB algorithms in scenario I. and lower bounds is known as the optimality gap.
It is said that the IAO approach converges to the proven optimality, when the optimality gap is less than the desired accuracy. This accuracy is tuned in this paper as 1e-10. The comparison of MCB and IAO methods is made, applying both methods on the same PC with an Intel(R) Xeon(R) CPU and 6 GB RAM, coded in Matlab/Yalmip environment and solved using Gurobi Optimizer 8.0. The restoration problem for all the test cases is considered as a multi-period optimization problem. The restorative period is assumed from 9:00 Am until 20:00 PM. The time step resolution is chosen to be 1 hour. We assume that the optimality threshold ( π πππ‘ ) and the computation time threshold ( π π‘πππ ) are set to 0.01 p.u. and 2 minutes, respectively. A. Scenario I: a small-scale off-outage area
The test network shown in Fig. 2 is a 11.4 kV balanced distribution network, which is based on a practical distribution system in Taiwan. It includes 2 substations, 11 feeders, 84 nodes, and 94 branches (incl. tie-branches). The detailed configuration data is given in [4]. Three dispatchable DGs on nodes 7, 39, and 80 have 2.8MW active and 3.0MVA apparent power capacities, while the DG on node 59 has 0.8MW active and 1.0MVA apparent power capacities. There exists also non-dispatchable DGs including a PV at node 28 and a Wind turbine at nodes 45. In scenario I, the restoration problem is solved for the test network shown in Fig. 2, where two faults occur at the same time on lines 30-31 (fault F1) and 43-44 (fault F2). These two faults are isolated by opening line switches {B5, 30-31} and {B6, 44-45}, respectively. The resulting off-outage area includes feeder E except node {30} and feeder F except nodes {43 and 44}. The restoration solution obtained from IAO and MCB approaches are reported in Table I. In order to deploy this solution on the network, first, the βReconfiguration Actionsβ must be implemented following the indicated order (I, II, etc.). Then, the βLoad Pickup Actionsβ are deployed throughout the subsequent instants of the restorative period according to the schedule given in Table I. Along with these results, Table I provides the optimal values of different objective terms and the computation time.
Fig. 5. The progress of the obtained solution using IAO and the MCB algorithms in scenario II.
The reliability objective term ( πΉ ππ ) is used to compare the quality of solutions provided by MCB and IAO methods . In this regard, the quality margin of a restoration solution is defined as the difference between its reliability objective value and the global optimal value of the reliability objective term. As reported in Table I, the solution provided by the proposed MCB approach is 13.75% far from the global optimal solution, provided by IAO method. The lower- and upper-bounds of the reliability objective term obtained using IAO and MCB approaches are plotted along their computation times in Fig. 4. In this regard, the lower and upper bounds of the MCB algorithm refer to the solutions provided, respectively, by the master and sub problem at each iteration. Whereas, for the IAO approach, each lower and upper bound correspond, respectively, to the best-bound and incumbent solutions found at a given iteration. It should be noted that in both methods only the upper bound solutions provide feasible solutions. As illustrated in Fig. 4, the best solution of the scenario I using the MCB algorithm is found at iteration 2 after 5.68 seconds. B. Scenario II: large-scale off-outage area in case of fault F3
In this scenario, the restoration problem is studied in case of fault F3 at substation 501 in the test network of Fig. 2. The off-outage area includes the whole feeders A, B, C, D, E, and F. In this case, the optimization problem P includes 529 binary variables including 23 reconfiguration variables ( y ) and 506 load pickup variables ( x ). As reported in Table I, it can be seen, the quality of the solution provided by the MCB method is 30.82% better than the one obtained with IAO method. In this scenario, IAO approach could not converge to the proven optimality. In this regard, the computation time that is given in Table I for IAO approach corresponds to the earliest time when IAO provides its best solution. In scenario II, the computation time threshold is met and we have to stop after the iteration number 2. However, for the sake of illustration goals, we let the iterations to continue until LB and UB solutions converge. The results are shown in Fig. 5. This figure shows the same type of information as illustrated in Fig. 4 but for test scenario II. As it can be seen, the quality of the best-found solution of the MCB method after 2 minutes is Since πΉ ππ has the largest weighting factor in the objective function, comparing the quality of the solution based on the overall objective value leads to the same result. Table II. Numerical results of test scenario III. S ce n a r i o S o l u ti on M e t hod R ec on f i gu r a ti on A c ti on s L o a d P i c kup A c ti on s πΉ π π ( p . u . ) πΉ π π€ ( m i n ) C o m pu t a ti on ti m e ( s ec ) III . a (f a u lt F ) M CB I. Open switches {6-7 ,11-12,25-26,27-28,28-29} and load breakers{1,3,5,19,24, 25} II.Close{T1,T4} III. Close load breakers {5,25}at 12:00 P.M. 75.5 334 61.12
III . b (f a u lt F ) M CB I. Open switches {31-32,33-40,41-42,42-43,59-60,63-64} and load breakers {32,33,34,35,52,53,54,55,58} II.Close{T1,T4} III. Close load breakers {35,58}at 12:00 P.M. 114.6 245.5 65.76
Fig. 6. The progress of the obtained solution using MCB algorithms in scenarios III.a and III.b. only 0.05% far from the quality of the final solution found at the convergence. The functionality of the IAO approach for solving scenario II is also illustrated in Fig. 5. As it can be seen, the optimality gap could not be reduced below 18% in IAO approach. C. Scenario III: large-scale off-outage area in case of faults F4 and F5
In order to further illustrate the performance of the proposed restoration algorithm, it is tested considering additional fault case studies. These are faults F4 and F5 represented in the test network of Fig. 3. This test network is based on a 11kV distribution network that is introduced in [39]. It includes 2 substations, 4 feeders, 70 nodes, and 76 branches (incl. tie-branches). One dispatchable DG with a capacity of 0.6 MW is installed on bus 68. The other type of DGs are PV-based DGs installed within LV networks at nodes 46, 47, and 61 with capacities of 0.6, 0.6 and 0.8 MW, respectively.
Table III. Checking network safety constraints for restoration solution obtained in different test cases S ce n a r i o Fault Min. voltage Margin (p.u.) Min. current margin (A) I F1 and F2 0.0236 p.u. at node 31 at time 11:00 A.M 9.716 A in line 84-11 at time 11:00 A.M II F3 0.0014 p.u. at node 29 at time 20:00 P.M. 12.82 A in line 85-47 at time 18:00 P.M. III.a F4 0.0013 p.u. at node 25 at time 14:00 P.M. 27.05 A in line 36-48 at time 13:00 P.M. III.b F5 0.0029 p.u. at node 31 at time 14:00 P.M. 19.04 A in line 67-68 at time 14:00 P.M.
In case of fault F4, feeders 1 and 2 will be in the off-outage area, whereas in case of fault F5, the off-outage area includes feeders 3 and 4. These two cases are presented, respectively, in scenarios III.a and III.b. The restoration problem contains 24 binary variables y and 186 binary variables x in case of scenario III.a; and 27 binary variables y and 234 binary variables x in case of scenario III.b. The IAO approach fails to present even a single feasible restoration solution for scenarios III.a and III.b. The progress of the MCB algorithm in solving the restoration problem for these scenarios is depicted in Fig. 6. The numerical results corresponding to the best solution found until the convergence of the MCB algorithm are given in Table II. It can be seen that it takes only 21 and 48 seconds for the MCB to find these best solutions in case of scenarios III.a and III.b, respectively. D. Discussion
The functionality of the proposed MCB with respect to the IAO method should be discussed separately for small scale-outage-areas such as in scenario I, and for large scale-outage-areas such as scenarios II and III. For smallβscale problems, the IAO method provides the optimal solution within a short time. As mentioned in section III.C, the proposed MCB algorithm does not necessarily provide the global optimal solution (as for IAO). However, as shown in scenario I, the quality of its solution is not far from the global optimal solution. The main advantage of the MCB method with respect to IAO method lies in large-scale optimization problems. Scenarios II and III illustrate this advantage. Actually, IAO approach could fail to converge to the optimal solution or even to provide a first feasible solution. As shown in scenario II, the best optimality gap obtained using IAO approach is significant and it means that the quality of the best feasible solution is poor. On the other hand, the MCB method found, within a very short time, a good-enough solution. There is no any unique and standard measure defining a good-enough restoration. However, it is obvious that the 0.05% quality margin that is obtained after the first iterations of the MCB algorithm in scenario II is acceptable and sufficient. In general, the IAO approach assigns the whole computational effort to finding the global optimal solution. If possible, this will be obtained for large-scale restoration problems after a long computation time. Whereas, the MCB method provides a good-enough solution at the first iterations thanks to the optimality cuts presented in section III.A. Then, through the subsequent
Table IV. Numerical restoration results obtained using the method of [24] in scenario II (fault F3 in the test network of Fig. 2). S o l u ti on M e t hod R ec on f i gu r a ti on A c ti on s L o a d P i c kup A c ti on s πΉ π π ( p . u . ) πΉ π π€ ( m i n ) C o m pu t a ti on ti m e ( s ec ) M e t hod o f [ ] _ t r y1 I. Open switches {5-6,12-13,38-39,44-45} and load breakers {12,14,17,19,26,27,28,29,31,32,34,35,36,37,38,41 ,46}, and close {T3,T7,T10,T11} II.Close {T1,T2,T4,T5,T8} III. Close load breaker 38 at 12:00 P.M. IV. Close load breakers {2,14,29} at 18:00 P.M. V. Close load breakers {35,36} at 19:00 P.M. 3.32e3 1.11e3 28.04 M e t hod o f [ ] _ t r y2 I. Open switch 33-34 and load breakers {2,3,4,11,12,14,16,17,19,20,21,26,27,28,29,31,32,33,34,35,36,37,38,40,41,42,44,46}, and close {T3,T7,T9,T11} II.Close {T2,T4, T8} III. Close load breakers {21,27,37} at 18:00 P.M. IV. Close load breakers {38,40} at 19:00 P.M. 6.32e3 1.23e3 18.82 iterations, it tries to improve the quality of the solution gradually. This characteristic is essential for the restoration problem, where an appropriate decision should be made in a very short time. In order to check the network safety constraints for each scenario, the obtained restoration solution is deployed on the model of the corresponding test network implemented in Matlab environment. The voltage and current profiles along the time are derived using power flow simulations in Matlab/MATPOWER toolbox. Table III gives the representative numerical results out of these profiles for each scenario. These results include the minimum nodal voltage magnitude and minimum line current margins over, respectively, all the nodes and lines of the networks and over all the time steps during the restorative period. These results show that the network safety constraints are all respected and therefore confirm the feasibility of the obtained solutions. Moreover, it can be seen that according to each restoration solution, the network is operated very close to its capacity envelop (especially in terms of the minimum voltage limit). This illustrates that within the safe region of the network operation, the most possible loads are restored for each scenario. E. Comparison with other mathematical programming methods
In this section, it is aimed to show the efficiency and superiority of the MCB with respect to the two mathematical programming methods proposed in [24] and [26]. In the first step, the MCB results in scenario II are compared with the results obtained from a mathematical formulation proposed in [24] for the restoration strategy. In [24], the electrical safety constraints are integrated into the restoration problem using DistFlow formulation. In this regard, the multi-period restoration problem is formulated in [24] as a mixed-integer linear programming model. According to this
Fig. 7. The voltage magnitude profile at different times steps during the restorative period (blue lines) and lower voltage limit (red dotted line) according to the solution obtained from the method of [24]_try1 in scenario II (nodes 43 and 44 are left without any supply according to the results of Table IV). formulation, the resulting topology of the network (i.e. line switching variables π¦ ππ ) could change at each time step. In order to make a fair comparison, we force the line switching variables in [24] to not change with time, as suggested in the proposed MCB methodology. With this modification, the formulation of [24] is implemented in Yalmip/Matlab and solved using Gurobi for the test case of scenario II. The obtained numerical results are reported in the first row of Table IV (Method of [24]_try1). The comparison of the πΉ ππ value in Table IV with πΉ ππ values reported in Table I for scenario II shows that the quality of the solution obtained using the method of [24] in try1 is better than the solution qualities of IAO and MCB approaches. However, since the DistFlow constraints are used in [24] instead of the AC power flow constraints, the obtained solution is infeasible regarding the minimum voltage limit. Fig. 7 confirms this infeasibility illustrating the results of a post power flow simulation for the obtained restoration solution. These results include the voltage magnitude profiles at different nodes and at different time steps during the restorative period. As it is shown, the lower voltage limit at some buses is violated at some time steps. The DistFlow approximation fails to guarantee the electrical safety constraints, since the network is operated close to its capacity envelop (i.e. current or/and voltage limits) during the restorative period. In case where voltage and/or line power constraints are violated, it is suggested in [24] to impose conservative limits on the DistFlow constraints. In this regard, we replace the original lower voltage limit (0.917 p.u.) to 0.970 p.u. This is the smallest value for the lower voltage limit in the DistFlow constraints that can guarantee the feasibility of the obtained solution in scenario II. Adding these conservative constraints, the restoration formulation of [24] is applied in try2 to solve the restoration problem in case of scenario II. The obtained results are shown in the second row of Table IV (Method of [24]_try2). This solution is feasible concerning all the electrical safety constraints. However, as it can be seen, the quality of the obtained solution (in terms of πΉ ππ ) is 52.74% lower than the Table V. Numerical restoration results in case of fault F6 in the test network of Fig. 2. S o l u ti on M e t hod S w it c h i ng A c ti on s πΉ π π ( p . u . ) πΉ π π€ ( m i n ) C o m pu t a ti on ti m e ( s ec ) M CB M e t hod I. Open switch 3-4 II.Close T2 4.4 60 4.04 M e t hod o f [ ] Step1 I. Open switch 6-7 II. Close {T1,T2} 1.65e3 150 1.14 Step2 Open load breakers {4,6} F6 DG DG OLTC2OLTC1 T1 B1 Feeder A B7 Feeder G B8 Feeder HT2 T12
B12 B13 B15B14
Critical and detachable loadsTie-switchClosed sectionalizing switchOpened sectionalizing switchClosed circuit breakerOpened circuit breaker
Fig. 8. The post-restoration configuration obtained in the first stage of the algorithm presented in [26]. solution quality of the proposed MCB method reported in Table I. This comparison clearly illustrates that it is not robust to simplify the AC power flow constraints in the restoration problem formulation with the corresponding DistFlow constraints. From the other hand, as the comparison of computation times related to the IAO and the method of [24] shows, the incorporation of AC power flow constraints increases the computation burden drastically. In order to make the restoration problem tractable in case of grids of realistic sizes while integrating the AC power flow constraints, the MCB method is proposed in this paper as a decomposition solution approach. In the next step, the proposed MCB methodology is compared with the restoration methodology presented in [26]. First of all, it should be mentioned that the methodology of [26] is particularly suitable for unbalanced networks and in case of extreme fault cases, where there is no access to the upper grid for the network restoration (referred as islanded network restoration). In this section, the comparison is conducted considering only assumptions made in this paper. It means that we focus on balanced distribution networks and we do not consider the islanded restoration of the distribution network. For this comparison, it is assumed that a fault occurs at the top point of feeder A (fault F6) in the test network of Fig. 2. This fault is isolated by opening the feeder breaker B1 and the switch on line 1-2. All the parameters of this test case are similar to the ones of scenario I and scenario II except that all the nodes are not equipped with load breakers. In this case study, it is assumed that only critical loads that are shown with β β β can be detached from their nodes. It means that among all the nodes in the off-outage area, only nodes 4 and 6 are equipped with load breakers. In order to make a fair comparison, we assume that load status variables ( πΌ π ) will not change with time, as proposed by the authors of [26]. With this modification, we apply the developed MCB formulation on this test case. The numerical results are reported in Table IV. According to this solution, since the switch on line 3-4 is opened, the nodes 2 and 3 are left without any supply. The restoration methodology of [26] is explained in section I. As mentioned, the authors of [26] propose to solve the restoration problem in two steps. In the first step, a heuristic approach is applied to find a suitable post-restoration topology for the network. This heuristic approach chooses a radial network configuration with the minimal diameter. The diameter of a tree is defined as the longest distance among all pairs of nodes in the network. The distance between a pair of nodes refers to the total impedance of lines on the shortest path between the two nodes. Applying this heuristic strategy to the test network of Fig. 2 in case of fault F6 results in the post-restoration configuration shown in Fig. 8. According to the restoration algorithm presented in [26], in the second stage, the status of load breakers and the outputs of sources are determined while fixing the network topology to the one obtained in the first stage. The resulting optimization problem is solved in [26] using a relaxed semi-definite programming methodology. But in this paper, we solve this same optimization problem using Gurobi solver, which adopts the branch-and-bound method to handle integrality constraints. The optimal solution is to open all the load breakers in the off-outage area (i.e. at nodes 4 and 6). This leads to a reliability objective term equal to 1.65e3. The comparison of πΉ ππ values given in Table IV shows that the quality of the solution obtained from the method of [26] is very far from the solution quality of the MCB method. It should be noted that if there were no load breaker in the off-outage area (no detachable loads), there would be no feasible solution for the optimization problem in the second stage of the algorithm presented in [26]. These numerical results clearly highlights the limits of [26] mentioned in section I. V. C ONCLUSION
The restoration is an NP-hard combinatorial optimization problem including three interdependent parts, namely, I) the reconfiguration problem, II) the load pickup problem, and III) the AC-Optimal Power Flow (AC-OPF) problem. This results in a huge and intractable problem especially considering a grid of realistic size in a multi-period problem. In order to relax the computation burden of the restoration problem, a two-stage decomposition approach is proposed in this paper named Modified Combinatorial Benders algorithm. In the outer level, the master problem solves a Mixed-Integer Linear optimization problem including the reconfiguration and the load pickup problems. The obtained configuration is fixed in the inner level and the load pickup variables are optimized subject to the AC-OPF constraints. The resulting sub problem is in the form of a Mixed-Integer Second-Order Cone Programming. This problem is broken down into several independent problems with smaller sizes. It makes the sub problem tractable in case of large-scale distribution networks. The solution of the sub problem is used to augment the feasibility or optimality cuts of the master problem. This algorithm is repeated through successive iterations until a solution with a desired level of optimality is obtained. The superiority of the proposed decomposition approach with respect to the integrated approach is illustrated with different scenarios on two test distribution networks. The results indicate the functionality of the Modified Combinatorial Benders method in providing, within a short time, a good-enough solution for large-scale restoration problems. The future major work will expand the proposed decomposition method in order to account for the uncertainties in the forecast of load demands and DG production. These uncertainties will be incorporated to the optimization problem resulting in a multi-period stochastic restoration problem. As another step forward, we plan to extend the proposed formulation to the unbalanced distribution network according to the approach that we sketched in Appendix A. VI. A CKNOWLEDGMENT
The authors gratefully acknowledge the financial support of the Qatar Environment and Energy Research Institute (QEERI). The authors are also warmly thankful to prof. Jean-Yves Le Boudec (EPFL, Lausanne β Switzerland) for his review regarding the mathematical formulation. VII. R EFERENCES [1] S. ΔurΔiΔ, C. S. Γzveren, L. Crowe, and P. K. L. Lo, βElectric power distribution network restoration: a survey of papers and a review of the restoration problem,β
Electr. Power Syst. Res. , vol. 35, no. 2, pp. 73β86, Nov. 1995, doi: 10.1016/0378-7796(95)00991-4. [2] H. D. de M. Braz and B. A. de Souza, βDistribution Network Reconfiguration Using Genetic Algorithms With Sequential Encoding: Subtractive and Additive Approaches,β
IEEE Trans. Power Syst. , vol. 26, no. 2, pp. 582β593, May 2011, doi: 10.1109/TPWRS.2010.2059051. [3] Y. Fu and H. Chiang, βToward Optimal Multiperiod Network Reconfiguration for Increasing the Hosting Capacity of Distribution Networks,β
IEEE Trans. Power Deliv. , vol. 33, no. 5, pp. 2294β2304, Oct. 2018, doi: 10.1109/TPWRD.2018.2801332. [4] C.-T. Su, C.-F. Chang, and J.-P. Chiou, βDistribution network reconfiguration for loss reduction by ant colony search algorithm,β
Electr. Power Syst. Res. , vol. 75, no. 2, pp. 190β199, Aug. 2005, doi: 10.1016/j.epsr.2005.03.002. [5] Y. Liu, R. Fan, and V. Terzija, βPower system restoration: a literature review from 2006 to 2016,β
J. Mod. Power Syst. Clean Energy , vol. 4, no. 3, pp. 332β341, Jul. 2016, doi: 10.1007/s40565-016-0219-2. [6] V. Gouin, M. A. Herault, and B. Raison, βStochastic integration of demand response and reconfiguration in distribution network expansion planning,β
Transm. Distrib. IET Gener. , vol. 12, no. 20, pp. 4536β4545, 2018, doi: 10.1049/iet-gtd.2018.5833. [7] I. A. Quadri, S. Bhowmick, and D. Joshi, βMulti-objective approach to maximise loadability of distribution networks by simultaneous reconfiguration and allocation of distributed energy resources,β
Transm. Distrib. IET Gener. , vol. 12, no. 21, pp. 5700β5712, 2018, doi: 10.1049/iet-gtd.2018.5618. [8] H. Sekhavatmanesh and R. Cherkaoui, βAnalytical Approach for Active Distribution Network Restoration Including Optimal Voltage Regulation,β
IEEE Trans. Power Syst. , pp. 1β1, 2018, doi: 10.1109/TPWRS.2018.2889241. [9] Lingwen Gan, Steven H. Low, βConvex relaxations and linear approximation for optimal power flow in multiphase radial networks,β presented at the Power Systems Computation Conference, Wroclaw, Poland, 2014, p. 9. [10] Emiliano DallβAnese, Georgios B. Giannakis, Bruce F. Wollenberg, βOptimization of unbalanced power distribution networks via semidefinite relaxation,β presented at the North American Power Symposium (NAPS), Champaign, IL, USA, 2012, p. 6. [11] Emiliano DallβAnese, Hao Zhu, Georgios B. Giannakis, βDistributed Optimal Power Flow for Smart Microgrids,β
IEEE Trans. Smart Grid , vol. 4, no. 3, pp. 1467β1475, Sep. 2013, doi: 10.1109/TSG.2013.2248175. [12] K. Chen, W. Wu, B. Zhang, and H. Sun, βSecurity evaluation for distribution power system using improved MIQCP based restoration strategy,β in
ISGT 2014 , 2014, pp. 1β6, doi: 10.1109/ISGT.2014.6816493. [13] N. C. Koutsoukis, D. O. Siagkas, P. S. Georgilakis, and N. D. Hatziargyriou, βOnline Reconfiguration of Active Distribution Networks for Maximum Integration of Distributed Generation,β
IEEE Trans. Autom. Sci. Eng. , vol. 14, no. 2, pp. 437β448, Apr. 2017, doi: 10.1109/TASE.2016.2628091. [14] J. A. Taylor and F. S. Hover, βConvex Models of Distribution System Reconfiguration,β
IEEE Trans. Power Syst. , vol. 27, no. 3, pp. 1407β1413, Aug. 2012, doi: 10.1109/TPWRS.2012.2184307. [15] R. A. Jabr, R. Singh, and B. C. Pal, βMinimum Loss Network Reconfiguration Using Mixed-Integer Convex Programming,β
IEEE Trans. Power Syst. , vol. 27, no. 2, pp. 1106β1115, May 2012, doi: 10.1109/TPWRS.2011.2180406. [16] M. Nick, R. Cherkaoui, J.-Y. L. Boudec, and M. Paolone, βAn Exact Convex Formulation of the Optimal Power Flow in Radial Distribution Networks Including Transverse Components,β
IEEE Trans. Autom. Control , vol. 63, no. 3, pp. 682β697, Mar. 2018, doi: 10.1109/TAC.2017.2722100. [17] H. Haghighat and B. Zeng, βDistribution System Reconfiguration Under Uncertain Load and Renewable Generation,β
IEEE Trans. Power Syst. , vol. 31, no. 4, pp. 2666β2675, Jul. 2016, doi: 10.1109/TPWRS.2015.2481508. [18] Z. Wang and J. Wang, βSelf-Healing Resilient Distribution Systems Based on Sectionalization Into Microgrids,β
IEEE Trans. Power Syst. , vol. 30, no. 6, pp. 3139β3149, Nov. 2015, doi: 10.1109/TPWRS.2015.2389753. [19] K. Chen, W. Wu, B. Zhang, and H. Sun, βRobust Restoration Decision-Making Model for Distribution Networks Based on Information Gap Decision Theory,β
IEEE Trans. Smart Grid , vol. 6, no. 2, pp. 587β597, Mar. 2015, doi: 10.1109/TSG.2014.2363100. [20] X. Chen, W. Wu, and B. Zhang, βRobust Restoration Method for Active Distribution Networks,β
IEEE Trans. Power Syst. , vol. 31, no. 5, pp. 4005β4015, Sep. 2016, doi: 10.1109/TPWRS.2015.2503426. [21] C. Lee, C. Liu, S. Mehrotra, and Z. Bie, βRobust Distribution Network Reconfiguration,β
IEEE Trans. Smart Grid , vol. 6, no. 2, pp. 836β842, Mar. 2015, doi: 10.1109/TSG.2014.2375160. [22] G. Li and M. Peng, βDynamic reconstruction of active distribution network with electric vehicles,β in , 2017, pp. 127β130, doi: 10.1109/ICCDS.2017.8120464. [23] A. Asrari, S. Lotfifard, and M. Ansari, βReconfiguration of Smart Distribution Systems With Time Varying Loads Using Parallel Computing,β
IEEE Trans. Smart Grid , vol. 7, no. 6, pp. 2713β2723, Nov. 2016, doi: 10.1109/TSG.2016.2530713. [24] B. Chen, C. Chen, J. Wang, and K. L. Butler-Purry, βMulti-Time Step Service Restoration for Advanced Distribution Systems and Microgrids,β
IEEE Trans. Smart Grid , vol. 9, no. 6, pp. 6793β6805, Nov. 2018, doi: 10.1109/TSG.2017.2723798. [25] H. Ahmadi and J. R. MartΓ, βDistribution System Optimization Based on a Linear Power-Flow Formulation,β
IEEE Trans. Power Deliv. , vol. 30, no. 1, pp. 25β33, Feb. 2015, doi: 10.1109/TPWRD.2014.2300854. [26] Y. Wang et al. , βCoordinating Multiple Sources for Service Restoration to Enhance Resilience of Distribution Systems,β
IEEE Trans. Smart Grid , vol. 10, no. 5, pp. 5781β5793, Sep. 2019, doi: 10.1109/TSG.2019.2891515. [27] J. Hooker,
Logic-Based Methods for Optimization: Combining Optimization and Constraint Satisfaction . John Wiley & Sons, 2011. [28] J. N. Hooker and H. Yan, βLogic circuit verification by Benders decompo-sition,β 1995. [29] G. W. Klau and P. Mutzel, βOptimal labeling of point features in rectangular labeling models,β
Math. Program. , vol. 94, no. 2, pp. 435β458, Jan. 2003, doi: 10.1007/s10107-002-0327-9. [30] G. Codato and M. Fischetti, βCombinatorial Bendersβ Cuts for Mixed-Integer Linear Programming,β
Oper. Res. , vol. 54, pp. 756β766, 2006, doi: 10.1287/opre.1060.0286. [31] A. P. Wierzbicki, βA mathematical basis for satisficing decision making,β
Math. Model. , vol. 3, no. 5, pp. 391β405, Jan. 1982, doi: 10.1016/0270-0255(82)90038-0. [32] Hossein Sekhavatmanesh, Rachid Cherkaoui, βA Multi-Step Reconfiguration Model for Active Distribution Network Restoration Integrating DG Start-Up Sequences,β accepted for publication in
IEEE Trans. Sustain. Energy , 2020, available online on Arxiv. [33] M. E. Baran and F. F. Wu, βOptimal capacitor placement on radial distribution systems,β
IEEE Trans. Power Deliv. , vol. 4, no. 1, pp. 725β734, Jan. 1989, doi: 10.1109/61.19265. [34] J. N. Hooker, βA Hybrid Method for Planning and Scheduling,β in
Principles and Practice of Constraint Programming β CP 2004 , 2004, pp. 305β316. [35] βStd, A. N. S. I. βC84. 1-2011.β American National Standard for Electric Power Systems and Equipment-Voltage Ratings (60 Hertz) (2011).β [36] E. Lopez, H. Opazo, L. Garcia, and P. Bastard, βOnline reconfiguration considering variability demand: applications to real networks,β
IEEE Trans. Power Syst. , vol. 19, no. 1, pp. 549β553, Feb. 2004, doi: 10.1109/TPWRS.2003.821447. [37] [Online], βHome Energy,β http://homerenergy.com , 15-Oct-2017. [Online]. Available: http://homerenergy.com. [38] [Online], βSolar Irradiation Data (SODA),β
Int. J. Electr. Power Energy Syst. , vol. 28, no. 5, pp. 331β338, Jun. 2006, doi: 10.1016/j.ijepes.2005.08.018. [40] Z. Liu and J. V. MilanoviΔ, βProbabilistic Estimation of Voltage Unbalance in MV Distribution Networks With Unbalanced Load,β
IEEE Trans. Power Deliv. , vol. 30, no. 2, pp. 693β703, Apr. 2015, doi: 10.1109/TPWRD.2014.2322391.
VIII. A PPENDICES A. Appendix I: Discussion on the Extension of the Model to Unbalanced Grids
As mentioned in section I, the proposed approach is specifically derived for balanced distribution networks. A potential way to extend this strategy to generic unbalanced grids is described in this appendix. It is assumed that the distribution network is operated under unbalanced but still radial conditions during the restorative period. Moreover, we neglect the impact of non-transposed or partially transposed (asymmetrical) lines. Since the lines in MV distribution networks are relatively short and have small impedances, the impact of line asymmetry on the voltage unbalance may be ignored as it is negligible with respect to the impact of unbalanced load and generation [40]. In this regard, all the electrical state variables (not the switching variables) including voltage, current, and power flow variables are decomposed using well-known sequence transformation method. As a result, the unbalanced grid is decomposed into three symmetrical and balanced three-phase circuits. To each of these circuits, we apply the relaxed-OPF formulation as given in section II. Regarding the transformation of the voltage/current limits from phase domain to the sequence domain, we apply the methodology given in [16]. In this regard, we make a conservative assumption. We assume that the negative and zero sequence terms of voltage and current magnitudes are binding to their standard/normal limits. Therefore, the voltage and current limits associated with the positive sequence terms are derived a priori . The constraints regarding the OPF formulation of three sequences are integrated into the Master- and Sub-problems. Once the optimization problem is solved, we transform the obtained values of electrical state variables from sequence domain back into the phase domain. B. Appendix II: Convex formulation of either-or constraints
As mentioned in section III.C, regarding the feasibility cut, at least one of two constraints (13) and (14) must hold. In order to integrate this either-or constraint into a convex optimization problem, binary variable π is introduced subject to the following constraints: β π¦ ππβπ(π,π¦ π) ) β€ (|π(π, π¦ (π) )| β 2) + π π (15) βπ (1 β π) + β π¦ ππβπ(π,π¦ π) ) β₯ 1 + β π¦ πππβπ(π,π¦ π) ) (16) where, π and π are two positive and sufficiently-large numbers. The auxiliary variable π determines which of the two constraints (13) and (14) must hold. According to (15), if π =0 , (13) is imposed and (14) is relaxed. When π = 1 , the situation is reversed. In both cases, one of the constraints is forced to be satisfied while the other constraint may also hold. C. Appendix III: Convex formulation of conditional constraints
Let βA implies Bβ denotes a conditional constraint, where A and B are logical expressions. This is logically equivalent to state that (
A and βΌ π΅ ) is false, where βΌ π΅ refers to the complement of B . The negation of ( A and βΌ π΅ ) is equivalent to ( βΌ A or π΅ ) . In section III.C, a conditional expression was given for the optimality constraint, which is stated again in the following: If (11) and (12) are satisfied, then β π π (1 β π₯ π ) πβπ(π,π¦ (π) ) β₯ β π (π¦ (π) ) must be satisfied. According to the above explanation, this is logically equivalent to the following expression: (13) or (14) or β π π (1 β π₯ π ) πβπ(π,π¦ (π) ) β₯ β π (π¦ (π) ) must hold As mentioned in section III.C, expressions (13) and (14) are the complements of (11) and (12), respectively. Therefore, the optimality cut can be translated into an either-or constraint and be reformulated according to the strategy given in section VIII.B. D. Appendix IV: Radiality Constraints
In this appendix, the radiality constraints formulated in (1.b)-(1.d) are briefly explained. The details and illustrations of these formulations can be found in [32]. The orientation of line ππ with respect to the available tie-switch (as virtual sources) is determined by continuous variables π½ ππ and π½ ππ . If the line is oriented from node π to node π , variable π½ ππ will be 1 and π½ ππ will be zero and if the line is oriented from node π to node π variable π½ ππ will be one and π½ ππ will be zero. In order to avoid possible isolated loops in cases where a DG exist in the off-outage area, (1.d) are added to the set of constraints. π is a large and positive coefficient. To each line ππ with switch in the out-of-service area, two continuous flow variables Ξ¨ ππ and Ξ¨ ππ are assigned. They are associated with the binary variables π½ ππ and π½ ππ , respectively. As formulated in the first two constraints of (1.d), for each line ππ with switch, at most one of the variables Ξ¨ ππ and Ξ¨ ππ gets a nonzero value depending on the line orientation that is identified with the variables π½ ππ and π½ ππ . The third constraint of (1.d) formulates the flow (in a virtual commodity) balance equation for each node ππ