1. Introduction

The future Smart Grid will be a flexible system that can quickly and reliably be reconfigured [5]. The ability of detecting and responding to system disturbance is an important characteristic of the Smart Grid network [17]. A large number of automated switching will be required to address the overload, voltage and phase imbalance issues caused by the presence of frequent plug-in electrical devices. An efficient distribution system reconfiguration strategy to quickly determine and optimize restoration plans can significantly improve the quality of the power supplying service. Therefore, a fully automated reconfiguration system is an important requirement of the Smart Grid which will be quick as well as reliable so that it should not compromise with the safety issue of the operating constraints of the whole network.

Network reconfiguration aims to modify the network structure of distribution feeders by changing the open/close status of the feeders (normally close) and tie (normally open) switches, attempting to reduce the general losses of the system by relieving the overloading of the network components, along with other objectives [23,4]. Optimal network reconfiguration for loss reduction is a highly complex combinatorial problem. This is caused by the large number of switching elements in a distribution network, and the nonlinear characteristics of the constraints that determine the electrical behaviour of the system. Furthermore, if we consider the other objectives like switching costs, etc, then we end up with a more difficult multiobjective problem. For this family of problems, several optimization techniques have been applied. In general, these optimization techniques can be separated into two large groups [6]: i) exact methods and ii) approximate methods. The exact methods attempt to find the global optimal solution by using algorithms similar to the branch and bound. However these methods can be applied only to simplified models of the electrical network with an approximated loss function. The approximate methods, on the other hand, use exact model of the network, and sometimes find a near-optimal solution. But there are no proof (or certificate) of optimality for these approximate methods; these are mostly heuristics.

The distribution system reconfiguration for loss reduction was first proposed in [23]. This method starts with a mesh network by considering all switches are closed. It then applies a recursive heuristic approach to open some switches in attempt to find a configuration with smaller operating loss. Baran and Wu [4] proposed another heuristic approach for loss reduction based on branch exchange. As expected, these search techniques, like any other heuristics, do not guarantee global optimization. Genetic algorithm (GA) often offers a better solution than the heuristic search methods, and takes less searching time [11]. In [11], the graph and matroid theories have been used to increase intelligence and efficiency of GA operators. Simulated annealing (SA) has also been used in [20,21] for network reconfiguration. SA is known to work well with nonconvex optimization problems, by helping the search algorithms to escape from local minima. The Ant Colony Search algorithm has also been proposed in [31] to find a near optimal solution for loss minimization. Many mixed integer programming algorithms, such as branch and bound, have also been utilized for distribution system reconfiguration [2,26,18]. In [26], an integer interior point method used a branch and bound technique to optimize the multi-objective system reconfiguration problem. In [18], the objectives of distribution system reconfiguration and optimal power flow are considered simultaneously and solve them jointly by using the Benders decomposition algorithm. Franco et al [12] developed an efficient mixed integer linear programming model for distribution system reconfiguration by considering the effect of distributed generation. The algorithm has been tested with both simulated and real radial distribution systems. The Enhanced Gravitational Search Algorithm (EGSA) has been applied for solving the reconfiguration problems with different objectives in [29,19]. The multi-objective reconfiguration framework developed in [29] considers the distribution line reliability along with two conventional reconfiguration objectives: operation cost and line loss minimization. The reliability factor is measured as a function of line failure rate, line length, average reconstruction time and reparation time. The effect of distributed generations on the transient stability of the distribution system has been considered in [19]. It has been shown in [33] that the critical clearing times (CCTs) associated with DG units can be used as a measure of the effect of DGs on the transient stability of system. Therefore, [19] incorporates CCT index of the system as a reconfiguration factor and solve the underlying multi-objective reconfiguration problem by using EGSA and fuzzy decision making method. The effect of variable loads on reconfiguration has been considered in [24] and the underlying problem is solved using an evolutionary GA. Stochastic method along with GA has been applied in [25] to reduce the power loss and operating cost of distribution system.

In this work, we model different objectives and constraints associated with a typical distribution system reconfiguration problem into a unified mathematical framework. We consider three objectives of distribution system reconfiguration: i) line loss minimization, ii) load balancing, and iii) switching cost minimization. At first we consider multiple possible paths to each bus from the substation. We represent them by a binary path matrix. We then develop an optimization algorithm that choose the best path from substation to each bus such that the objectives are minimized. One of the principal problems arises when we incorporate the so called radiality constraint within the mathematical formulation. This radiality constraint ensures that the network topology is radial. The demand to ensure a radial network topology is hard to incorporate into the mathematical model in an efficient manner. This limits the applicability of the existing methods significantly. We handle the constraint by using some binary variables. The resulting optimization problem has a quadratic form which we solve efficiently by using mixed integer programming.

2. Problem statement

2.1. Basic power flow equations

Let us consider the radial distribution network in Figure 1. The line impedance of a branch l is represented by ${z_l} = {r_l} + {\text{i}}{x_l}$, and the load at bus i is ${S_{Ni}} = {P_{Ni}} + {\text{i}}{Q_{Ni}}$. Power flow in the radial network can be described by a set of DistFlow branch equations [4]:

$P_{i+1}=P_i-r_i\frac{P_i^2+Q_i^2}{V_i^2}-P_{Ni+1},$

(1)

$Q_{i+1}=Q_i-x_i\frac{P_i^2+Q_i^2}{V_i^2}-Q_{Ni+1},$

(2)

$V_{i+1}^2=V_i^2-2(r_iP_i+x_iQ_i)+(r_i^2+x_i^2)\frac{P_i^2+Q_i^2}{V_i^2},$

(3)

where V_i is the voltage of bus i. The quadratic terms in the equations (1)-(3) represent the line loss which are smaller compared to line power P_i and Q_i. Therefore, by dropping the second order terms we can produce simplified DistFlow equations [4]:

$P_{i+1}=\sum_{j=i+2}^nP_{Nj},$

(4)

$Q_{i+1}=\sum_{j=i+2}^nQ_{Nj},$

(5)

$V_{i+1}^2=V_i^2-2(r_iP_i+x_iQ_i).$

(6)

Also the branch power loss can be approximated as

$L_i=r_i\frac{P_i^2+Q_i^2}{V_i^2}=\bar{r}_i(P_i^2+Q_i^2). \ \mathrm{p.u.}$

(7)

where $\bar{r}_i=\frac{r_i}{V_i^2}$.

2.2. Objective terms

Often a distribution system is reconfigured to minimize the total branch power loss of the system:

$L=\sum_{i=0}^nL_i$

(8)

To avoid equipment overload, and ensure the power supply quality, another general objective considered for system reconfiguration is load balancing [4,7]. Let S _i be the complex power at the sending end of branch i, and $S_i^{\max}$ be the KVA capacity of the branch. Then one defines the load balancing index of branch i as [4,31]:

$G_i=\left(\frac{S_i}{S_i^{\max}}\right)^2=\frac{P_i^2+Q_i^2}{(S_i^{\max})^2}.$

(9)

Clearly, G_i is a quantitive measure of the extent to which the branch i is loaded. Then the objective of load balancing of the whole system would be to minimize the criterion [7] :

$\max_{i=0, \cdots n}G_i.$

(10)

Note that minimizing the criterion(10) results into a minimax problem.

Another target of distribution system reconfiguration is to restore all consumer service as quickly as possible. Many factors affect the system restoration time. In many cases, the restoration time can be improved by minimizing switching operating cost [8,28]. Let s_i be a binary variable, whose value will be 1 if the state of a switch at branch i is changed in reconfiguration. For example, if the switch at branch i was normally closed, however, the system reconfiguration demands to open it, then s_i = 1. Otherwise, if the system reconfiguration demands the switch to remain close then s_i = 0. Let c_i be the cost associated with the state changing of s_i. Then the objective of switching cost minimization can be expressed as to choose s such that

$\sum_{i=0, \cdots, n} c_i s_i$

(11)

is minimized.

Along with these objectives we need to maintain few constraints [31]:

1. Voltage magnitude of each bus must lie within its permissible range.

2. Current magnitude of each branch must lie within its permissible range.

3. Total load of a particular network can not exceed the capacity of the corresponding power source.

4. Distribution networks should be composed of radial structure operation.

3. Problem formulation

3.1. Path matrix

Let us consider a distribution network with n branches (including tie branches), $\hat{m}$ substations and m distribution buses. Assume that every branch can be closed or opened by a switch placed on the branch. A path is defined as a directed sequence of connected branches that connects two buses. For example, in Figure 2, Bus 0 is a substation and Branch 8 is a tie branch. Hence, we have n = 8; $\hat{m}$ = 1; m = 7. A possible path from Bus 0 to Bus 3 can be Branch1→Branch2→Branch3. Given a bus l, let $v^{(l)}\in \mathbb{R}^n$ be a binary vector such that the value of $v^{(l)}_i$ represents the status of the switch on the branch i of the network. If $v^{(l)}_i$ = 1 then the switch on the branch i is closed. Then we use the binary vector $v^{(l)}$ to represent a particular switching combination of all branches so that the bus l can be reached from a substation. In other words, a path will be created from the substation to bus l by “selecting the switch states” in $v^{(l)}$. There may be multiple paths between the substation/s and bus l. Suppose, there are ql possible paths between the substation/s and bus l. Then the q_l switching combinations are indicated by $\{v^{(l, j)}\}_{j=1}^{q_l}$. To illustrate with an example let us consider Figure 2 where Bus 0 is the substation. There are two possible paths between Bus 0 and Bus 2, hence q₂ = 2 and

$\begin{align} v^{(2, 1)}=[1 \ 1 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0]^T\nonumber\\ v^{(2, 2)}=[1 \ 0 \ 1 \ 0 \ 1 \ 1 \ 0 \ 1]^T\nonumber \end{align}$

Let us construct m submatrices B⁽¹⁾, B⁽²⁾, ... B^(m) such that $B^{(l)}=[v^{(l, 1)} \ v^{(l, 2)} \cdots v^{(l, q_l)}]$. Then every column of B^(l) indicates a possible path between the substation and bus l. Now to reach all buses from the substation in the distribution network we need to select exactly one column from each submatrix $\{B^{(l)}\}_{l=1}^m$. Note that, by selecting more than one columns from a submatrix, we create one or more self loop/s in the system, thereby violating the required radial topology.

3.2. Radial constraint

Different methods can be used to avoid loop in the network topology. For example, the procedure used in [32,3,27] construct an initial network by closing all tie and branch switches of the distribution system. Then they count the number of switches in every potential loops in the network. The main principle of avoiding loop in the final distribution network is to ensure that there should be at least one open switch in each potential loop of the network. The procedure needs to know the status of each switch to ensure the radial constraint. In this work, we shall develop a methodology that can maintain radial constraint of the network by using the path matrix B.

To enforce the the radial topology mathematically, let us construct

$B=[B^{(1)}, B^{(2)}, \cdots B^{(m)}],$

where $B\in \mathbb{R}^{n \times M}$ and $M=q_1+q_2+\cdots q_m$. Let $z\in \mathbb{R}^M$ be a binary vector which is acting as a column selector of B. If $z_i=1$ then the i-th column of B is selected. For example, in Figure 2 we can construct,

$B = \left[ {\underbrace {\begin{array}{*{20}{c}} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{array}}_{{B^{\left( 1 \right)}}}\underbrace {\begin{array}{*{20}{c}} {11} \\ {10} \\ {01} \\ {00} \\ {01} \\ {01} \\ {00} \\ {01} \end{array}}_{{B^{\left( 2 \right)}}}\underbrace {\begin{array}{*{20}{c}} {11} \\ {10} \\ {10} \\ {00} \\ {01} \\ {01} \\ {00} \\ {01} \end{array}}_{{B^{\left( 3 \right)}}}\underbrace {\begin{array}{*{20}{c}} {11} \\ {10} \\ {10} \\ {11} \\ {01} \\ {01} \\ {00} \\ {01} \end{array}}_{{B^{\left( 4 \right)}}}\underbrace {\begin{array}{*{20}{c}} {11} \\ {10} \\ {10} \\ {00} \\ {01} \\ {10} \\ {00} \\ {10} \end{array}}_{{B^{\left( 5 \right)}}}\underbrace {\begin{array}{*{20}{c}} {11} \\ {10} \\ {10} \\ {00} \\ {01} \\ {01} \\ {00} \\ {10} \end{array}}_{{B^{\left( 6 \right)}}}\underbrace {\begin{array}{*{20}{c}} {11} \\ {10} \\ {10} \\ {00} \\ {01} \\ {01} \\ {11} \\ {10} \end{array}}_{{B^{\left( 7 \right)}}}} \right]$

(12)

Since the test system in Figure 2 has m = 7 buses, the matrix B is concatenation of 7 submatrices $B=[B^{(1)}, B^{(2)}, \cdots B^{(7)}]$. Note that the matrix B has total M = 13 column and every columns represents a possible path between the substation node to a particular bus. Therefore, the vector $z$ has 13 components. Now if $z_3=1$ then the 3rd column of B will be selected and hence we choose the path Branch1→Branch5→Branch6→Branch8→Branch3 to reach Bus 2.

Let us partition $z$ as

$z = [z^{(1)} \\ z^{(2)} \\ \cdots \\ z^{(m)}]^T$

(13)

where the vector $z^{(l)}$ is of length q_l. To avoid creating self loops and ensure every bus is reachable we need to impose the following constraints on $z$:

$\sum_{i=1}^{q_l} z_i^{(l)} =1; \ \ \mathrm{for} \ l=1, \cdots m.$

(14)

However, this constraint is not sufficient to rule out the possibility of generating a loop by two paths to two different buses. For example, in (12), $z_1^{(2)}=0, z_2^{(2)}=1, z_1^{(3)}=1$ and $z_2^{(3)}=0$. Having $z_2^{(2)}$ will select the 3rd column of B, and it allows to reach Bus 2 via the path Branch1→Branch5→Branch6→Branch8→Branch3, and $z_1^{(3)}=1$ will select the 4th column of B, and it allows to reach Bus 3 via the path Branch1→Branch2→Branch3. The constraint in (14) is satisfied. However, those paths together create a loop Branch1→Branch5→Branch6→Branch8→Branch3→Branch2→Branch1. To avoid such situations, we need addition constraints.

To deal with the above issue we introduce the notation of a “subpath”. We say the path associated with a column $v^{(l, j)}$ is a subpath of the path associated with $v^{(i, k)}$ if

${v^{(i,k)}}{ \odot ^{(l,j)}} = {v^{(l,j)}},$

(15)

where $\odot$ denotes the component wise multiplication operator between two vectors. A close inspection to (15) reveals that $v^{(l, j)}$ will be a subpath of $v^{(i, k)}$ if

1. Both paths start from the same node.

2. Buses visited by path $v^{(l, j)}$ are also visited by $v^{(i, k)}$. Note that $v^{(i, k)}$ may visit additional buses which are not visited by $v^{(l, j)}$.

3. If the status of a particular branch is ON in path $v^{(l, j)}$, then the status of the branch is also ON in $v^{(i, k)}$.

4. If $i\not=j$, then by construction, for a fixed l, $v^{(l, j)}$ cannot be a subpath of $v^{(l, j)}$.

If (15) holds, we express this as

$v^{(i, k)} \supset v^{(l, j)}.$

For these rules it is clear that, if the path $v^{(l, j)}$ is a subpath of $v^{(i, k)}$, and the path $v^{(i, k)}$ is activated, then $v^{(l, j)}$ is activated as well.

To illustrate the idea with an example, consider again Figure 2 and the associated B matrix in (12). Suppose we make $z_1^{(3)} = 1$ i.e, we choose $v^{(3, 1)}$. This will select the 4th column of B, creating a path Branch1→Branch2→Branch3 to reach Bus 3. Now to reach Bus 2 we can choose either $v^{(2, 1)}$ (2nd column of B) or $v^{(2, 2)}$ (3rd column of B). But it is readily verified that $v^{(2, 1)} \subset v^{(3, 1)}$, and $v^{(3, 1)}$ is already selected. Hence $v^{(2, 1)}$ will be selected automatically, so that we reach Bus 2 via Branch1→Branch2. In our optimization algorithm we can easily account for the above constraint. When a particular path is selected, all its subpaths must also be selected. Mathematically, we have the constraint

${v^{(i,k)}} \supset {v^{(l,j)}} \qquad \Rightarrow z_k^{(i)} ≤ z_j^{(l)}.$

(16)

3.3. Branch power loss

For every bus l, construct a diagonal matrix

$U^{(l)}=S_{Nl}I_{q_l\times q_l}.$

Then form

$U = \left[ {\begin{array}{*{20}{c}} {{U^{\left( 1 \right)}}}&{{0_{{q_1} \times {q_2}}}}& \cdots &{{0_{{q_1} \times {q_m}}}} \\ {{0_{{q_2} \times {q_1}}}}&{{U^{\left( 2 \right)}}}& \cdots &{{0_{{q_2} \times {q_m}}}} \\ \vdots & \vdots & \cdots & \vdots \\ {{0_{{q_m} \times {q_1}}}}&{{0_{{q_m} \times {q_2}}}}& \cdots &{{U^{\left( m \right)}}} \end{array}} \right]$

(17)

Let us define

$H=BU.$

Suppose we have chosen a $z$ satisfying the constraints (14) and (16) to reach all buses of the n branch m bus radial network. After a few steps of algebra using the simplified DistFlow equations in (4)-(6), one can verify that the l th component of the product term $Hz$ is equal to power flow via the l th branch. For the example in Figure 2 if we take $z=[1 \ 1 \ 0 \ 1 \ 0 \ 1 \ 0 \ 0 \ 1 \ 0 \ 1 0 \ 1]^T$, then:

$Hz = \left[ {\begin{array}{*{20}{c}} {\sum\nolimits_{i = 1}^7 {{S_{Ni}}} } \\ {\sum\nolimits_{i = 2}^4 {{S_{Ni}}} } \\ {\sum\nolimits_{i = 3}^4 {{S_{Ni}}} } \\ \vdots \\ 0 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{P_1} + {\text{i}}{Q_1}} \\ {{P_2} + {\text{i}}{Q_2}} \\ {{P_3} + {\text{i}}{Q_3}} \\ {{P_4} + {\text{i}}{Q_4}} \\ \vdots \\ {{P_8} + {\text{i}}{Q_8}} \end{array}} \right].$

(18)

Let R be a n × n diagonal matrix where the diagonal components are the branch resistance, i.e.,

$R_{i, i}=\bar{r}_i; \ \mathrm{for} \ i=1, \cdots, n.$

If the constraints (14) and (16) hold, then using (7) and (8), the total line loss is given by

$(Hz)^* R (Hz) = z^*Wz,$

(19)

where

$W = H^*RH$

is a known complex-valued Hermitian non-negative definite matrix; it has the form W_r + iW_i, where W_r is real valued symmetric matrix, and W_i is a real valued skew symmetric matrix. Because z is real valued, ${z^*}{W_{\text{i}}}z = 0$, implying

${z^*}Wz = {z^*}{W_{\text{r}}}z.$

Now the first objective in (8) can be written as

$\mathop {min}\limits_z {z^*}{W_{\text{r}}}z$

(20)

subject to $z_i=\{0, 1\}$, and constraints (14), (16). Note that optimizing (20) is a nonlinear problem because R is a function of V. However, if we fix the value of V, then the problem can be solved by using a standard quadratic binary programming solver [1]. In Table 3 we have developed a two step iterative algorithm to solve the optimization problem in (20).

3.4. Load balancing

For load balancing, the problem of minimizing the criterion (10) can be formulated as

$\begin{gathered} \min y \\ {\text{subject}}\;{\text{to}},\;{G_i} ≤ y;\;\;\forall \;i \\ y > 0 \\ \end{gathered}$

(21)

with the constraints (14) and (16). Recall that the l th component of $Hz$ is equal to the power flow via the l th branch, see the discussion around (18). Hence, according to (9),

$G_i=(H(i, :)z)^T\frac{1}{(S_i^{\max})^2}(H(i, :)z).$

(22)

Then (21) can be written as

$\mathop {min}\limits_{z,y} y$

(23)

$\begin{gathered} {\text{subject}}\;{\text{to}},\;{(H(i,:)z)^T}\frac{1}{{{{(S_i^{\max })}^2}}}(H(i,:)z) ≤ y;\;\;\forall \;i \\ y > 0 \\ \end{gathered}$

(24)

3.5. Switching cost

Recall from (11) that the switching cost depends on a binary vector s, where s_i = 1 when the state of a switch at branch i is changed in reconfiguration. So far we have formulated all optimization problems in terms of the variable $z$. To extend our framework to incorporate the switching costs, we first define a linear map from $z$ to s. First we express the binary vector of branch switch status in terms of $z$. In particular,

1. every branch switch of the network is mapped by at least one element of $z$ i.e., $z_i$,

2. every $z_i$ must be mapped to exactly one branch switch.

As described in Section 3.1, every row of the matrix B corresponds to a particular branch switch. Let us define the length $\mathcal{L} \{ v^{(l, j)} \}$ of the path generated by a column $v^{(l, j)}$ as the number of branches in that path. In particular, if 1 denotes the M dimensional vector of all ones, then

$\mathcal{L} \{ v^{(l, j)} \} = {\bf 1}^T v^{(l, j)}.$

Having $z^{(l)}_j = 1$ means that the path associated to $v^{(l, j)}$ is created by closing all the branch switches on the path. If $\mathcal{L}\{ v^{(l, j)} \}>1$, there is always a bus l₁ and an index j₁ such that $v^{(l_1, j_1)} \subset v^{(l, j)}$, and

$\mathcal{L}\{ v^{ ( l_1, j_{1} ) } \}=\mathcal{L}\{ v^{(l, j)} \} - 1,$

(25)

meaning $v^{(l, j)}$ has only one extra branch than $\{ v^{(l_1, j_1)} \}$. Again, the constraint in (16) implies that when $z^{(l)}_j = 1$, i.e. $v^{(l, j)}$ is selected, then all its subpaths must also be selected. Consequently, if $z^{(l)}_j = 1$, then $z^{(l_1)}_{j_1} = 1$. Having $z^{(l_1)}_{j_1} = 1$ requires to close all switches on $v^{(l, j)}$, except the only switch on the branch which is not on $v^{(l_1, j_1)}$. Similarly, if $\mathcal{L} \{ z^{(l_1, j_1)} \} > 1$, then we can find a bus k₂ and an index j₂ such that $v^{(l_2, j_2)} \subset v^{(l_1, j_1)}$, and

$\mathcal{L}\{ v^{ ( l_2, j_{2} ) } \}=\mathcal{L}\{ v^{(l_1, j_1)} \} - 1,$

Again, having $z^{l_2}_{j_2} = 1$ requires to close all switches on $v^{(l_1, j_1)}$, except the only switch on the branch which is not on $v^{(l_2, j_2)}$. The same process continues until the length of the path is one. Thus, in process of creating the path $v^{ ( l, j ) }$, the variable $z^{(l)}_j$ need not be made responsible to close all branch switches. Instead, it is suficient that having $z^{(l)}_j = 1$ is equivalent of closing the switch in the branch which is in $v^{(l, j)}$, but is not present in any of its subpaths. For example, in (12), $v^{(3, 1)}$ (4th column of B) is associated to the path Branch1→Branch2→Branch3. Hence, when $z^{(3)}_1=1$ we select $v^{(3, 1)}$, and we need to close switches on all of these three branches. Again, $v^{(1, 1)} \subset v^{(2, 1)} \subset v^{(3, 1)}$. Now Branch3 is the only branch in $v^{(3, 1)}$, which does not belong to any of its subpaths. Similarly, Branch2 is the only branch in $v^{(2, 1)}$, which does not belong to any of its subpaths, while the only branch in $v^{(1, 1)}$ is Branch1. According to our proposal, making $z^{(3)}_1 = 1$ closes the switch in Branch3. But, by (16), $z^{(3)}_1 = 1$ implies $z^{(2)}_1 = 1$, and thus the switch on Branch2 will have to be closed by $z^{(2)}_1 = 1$. As this inductive process continues, eventually, all the switches in the path Branch1→Branch2→Branch3 are closed.

Based on the above idea, we develop a generic algorithm in Table 2 to generate a mapping matrix D such that the product $D z$ gives the switch status in the network. In particular, $D z$ is a binary vector, where, if the j th component of $D z$ is one, then the switch on the j th branch needs to be closed, otherwise the switch should be opened. For the radial network in Figure 2, the mapping matrix D, and switch status of different branches for $z=[1 \ 1 \ 0 \ 1 \ 0 \ 1 \ 0 \ 0 \ 1 \ 0 \ 1 0 \ 1]^T$ is

$\underbrace {\left[ {\begin{array}{*{20}{c}} 1 \\ 1 \\ 1 \\ 1 \\ \vdots \\ 0 \end{array}} \right]}_{S witch} = \underbrace {\left[ {\begin{array}{*{20}{c}} 1&{0 0}&{0 0}&{0 0}&{0 0}&{0 0}&{0 0} \\ 0&{1 0}&{0 0}&{0 0}&{0 0}&{0 0}&{0 0} \\ 0&{0 1}&{1 0}&{0 0}&{0 0}&{0 0}&{0 0} \\ 0&{0 0}&{0 0}&{1 1}&{0 0}&{0 0}&{0 0} \\ 0&{0 0}&{0 0}&{0 0}&{0 1}&{0 0}&{0 0} \\ 0&{0 0}&{0 0}&{0 0}&{1 0}&{0 1}&{0 0} \\ 0&{0 0}&{0 0}&{0 0}&{0 0}&{0 0}&{1 1} \\ 0&{0 0}&{0 1}&{0 0}&{0 0}&{1 0}&{0 0} \end{array}} \right]}_Dz$

(26)

Let the state of the switches before reconfiguration be represented by a binary vector d, where d_i = 1 denotes the switch in the branch i is closed. Now construct a vector e such that

$e_i=\left\{\begin{array}{cc} -1, & d_i=1 \\ 1, & \mathrm{otherwise}.\end{array} \right.$

(27)

Now the change of switching state after distribution system reconfiguration can be represented as

$s=d+EDz$

(28)

where $E=\mathrm{diag}\{e\}$ and s as in (11).

3.6. Current constraint

As described in Section 2.2, current magnitude of each branch of the final distribution network must lie within its permissible range. Let $I_l^{\min}, I_l^{\max}$ are the magnitudes of minimum and maximum permissible current in the branch l. Then the current constraint can be written as

$(I_l^{\min})^2\le(H(l, :)z)^T\frac{1}{V_l^2}(H(l, :)z)\le (I_l^{\max})^2;\forall l$

(29)

3.7. Final optimization problem

By combining all objectives the final optimization problem can be represented as

$\min_{z, y} \ \left\{\tau_1(z^*W_rz)+\tau_2 y +\tau_3(c ED z)\right\}$

(30)

subject to $z_i\in \{0, 1\}$ and the constraints in (14), (16), (24), (29). Where $\tau_i\ge0$ is related to the weighting of different objectives. We can use a mixed integer quadratic programming (MIQP) solver, like GUROBI optimization solver [1], to optimize (30). The final algorithm is given in Table 3. The algorithm starts with assuming V_i = 1 p.u. After solving (30) in Step 1, the algorithm construct a distribution network using resultant z. The voltage and current profile of new network is calculated in Step 3. If the voltages remain their permissible range, then the algorithm update $V_i^{(l)}$ with the new V_i and go to Step 5. If the voltage constraints are not satisfied then the current solution is not feasible, and hence the algorithm adds a constraint with (30) to ensure that the current solution will not come again (Step 4 b)). Step 5 checks whether two successive solutions are same. If two solutions are same then we get the optimum solution, and the algorithm terminates. Otherwise the algorithm returns to Step 1.

Note that the size of the path matrix B in (12) increases with increasing the number of tie switches in the distribution network. The interesting fact is that the radial configuration of distribution system corresponds to a “spanning tree” of a graph representing the network topology [4]. Our problem is to seek a minimal spanning tree that satisfies the objective and constraints in (30) [4]. The experiment result shows that the optimization algorithm always try to select a smaller path (i.e., a path with smaller number of branches) to reach from a substation to a particular bus. Hence, while constructing B^(l) in (12) for large network, we can choose N number of smallest possible paths from the substation/s to bus l. In this work, we utilize Yen’s k-Shortest path algorithm [34] to find first N shortest paths between two nodes. In our simulations, we set N = 10.

4. Experimental results

The performance of the proposed method is evaluated by reconfiguring three distribution test systems: 32-bus [4], 70-bus [9] and 135-bus [16]. The resulting mixed integer quadratic programs in (30) were solved by using the solver GUROBI [1]. The objective function was modeled by using MATLAB programming. Load flow solutions were computed using MATPOWER [35]. For load balancing, we define a balancing index term $\mathbb{G}$ such that

$\mathbb{G}=\max_i \frac{P_i^2+Q_i^2}{(S_i^{\max})^2}.$

(31)

4.1. 32-Bus System

The 32-bus distribution system is shown in Figure 3. The system consists of one source transformer, 32-bus-bars, and 5 tie switches. The detail data of line impedance and load demand can be found in [4]. The total active and reactive power for the whole system loads are 5048.26 kW and 2547.32 kVAR, respectively. The current carrying capacity of branch No.1 - 9 are 400A, and the other remaining branches including the tie lines are 200A [10].In the following, a line switch between bus j and k is denoted by (j, k).

Most of the conventional algorithms does not consider optimizing multiobjectives (loss minimization, load balancing and switching) at the same time. Hence, for a fair comparison of the proposed algorithm with other algorithms, at first, we only consider loss minimization. We use the optimization in (30) where we set $\tau_1=1$ and $\tau_2=0$, $\tau_3=0$. The initial losses of the test system was 202:68 kW. Table 4 shows a comparison between various algorithms. As can be seen, the methods of Goswami [15], Gomes [13], McDermott [22] find the optimum configuration. The Shirmohammadi [30] method also find a near-optimum solution. In the case, we found that most of the algorithms have obtained near-optimum solution since the system is small. Table 4 also shows the computation times required by different algorithms for solving the optimization problem. The simulations are performed in MATLAB using an Intel Core 2 Duo 2.66 GHz processor with 4-GB of memory.

We then overload some buses of the test system and investigate the performance of different algorithms for loss minimization. The output will demonstrate the robustness of different algorithms. We consider the overloading proposed in [14]. The loads at bus 9 and bus 13 were 60 kW 20 kVAR and 120 kW 80 kVAR respectively. The load of both buses were changed to 420 kW 200 kVAR. The initial line losses become 299.94 kW. Table 5 shows a comparison of system reconfiguration performance of different algorithms. As demonstrated in [14], the optimum results obtained by using the exhaustive evaluation algorithm suggests to open (8, 9), (13, 14), (27, 28), (31, 32), (7, 20). Hence, the proposed algorithm gets the optimum solution. It saves 32:9% the total line loss. Another interesting observation is that although McDermott [22] method found optimum solution in Table 4, however, it shows worst result in Table 5. Table 6 shows the currents of two main feeders before and after system reconfiguration.

We include the switching cost to the optimization in (30). As demonstrated in [8], we set the opening cost of a feeder switch is 1 and closing cost of a tie switch is 0.5. In particular, we set

$c_i=\left\{\begin{array}{cc} 1, & \mathrm{If \ branch \ i \ is \ a \ line \ feeder} \\ 0.5, & \mathrm{if \ branch \ i \ is \ a \ tie \ feeder}.\end{array} \right.$

(32)

We consider the actual 32-bus system [4] to illustrate the effect of multiobjective optimization on system reconfiguration. We set different values of $\tau_i$ (30) to check the performance in Table 7. As can be seen, the switching cost increases with increasing the loss savings. Another interesting issue is that if we change the value of $\tau_3$ from 0.004 to 0.001, the switching cost jumps from 3 to 7.5, whereas the loss savings increases only 2.4%. This may be an important consideration issue while reconfiguring distribution system. The effect of load balancing can be visualized effectively when we consider a distribution system where power is supplied from multiple feeders. In the following, we consider some multifeeder distribution systems.

4.2. 70-bus system

The 70-bus test system is shown in Figure 4. It is an 11 kV distribution system that has two substations, 70 buses, 68 feeder lines and 8 tie lines. System data can be found in [9]. The base voltage and power of the system are 11 kV and 10 MVA respectively. Before network reconfiguration, total loss of the system was 341:43 kW. Table 9 shows the reconfiguration results for loss minimization by different algorithms. In the optimization of (30), we set $\tau_1=1$ and $\tau_2=0$, $\tau_3=0$. The proposed method recommends to open the switches (9, 15), (21, 27), (28, 29), (43, 38), (40, 44), (49, 50), (62, 65), (67, 15). As can be seen the proposed algorithm can performs better than [9]. We then consider multiobjective optimization in Table 8. As before, we set the values of ci as in (32). There are some interesting observations. Let us consider the case $\tau_1=1$, $\tau_3=0.01$, where we do not consider load balancing i.e., $\tau_2=0$. The switching cost is 1.5, however the load balancing index is high =0.5476. Again, when we take $\tau_2=1$, the line loss remain almost same, the value of load balancing index decreases significantly ($\mathbb{G}=0.3646$), however switching cost increases to 3. Hence, one case calibrate the value of $\tau_i$ to reach a desire target.

4.3. 135-Bus System

To illustrate the performance of the proposed framework in a relatively large network we consider an 13.8 kV distribution system that has eight substations, 143 buses, 135 branch lines and 21 tie lines [16]. System data can be found in “http://venus.ece.ndsu.nodak.edu/～kavasseri/reds.html”. The base voltage and power of the system are 13.8 kV and 10 MVA respectively. The reconfiguration results for different settings are shown in Table 10.

5. Conclusion

The paper presents a mathematical framework for distribution system reconfiguration. The optimization formulation is in quadratic mixed integer form, and we solve them efficiently by using a solver called GUROBI. In experimental section we demonstrate different aspects of system reconfigurations. The proposed approach considers some practical objectives associated with distribution system reconfiguration problem and applied it efficiently to reconfigure some realistic test systems. Furthermore, the algorithm is also capable of reconfiguring relatively large scale distribution system such as 135-bus test system. Therefore, the proposed methodology can be applied in practical distribution system reconfiguration.

Although, the proposed method considers many important reconfiguration factors like brach power loss minimization, load balancing, switching cost, there exists some others factors that can improve the reconfiguration performance of real distribution system. For example considering the effect of adding/removing distributed generations (DGs) on distribution system is a factor that can improve the reconfiguration performance [19,12]. The proposed mathematical model is general enough to incorporate the effect of changing DG connections on the reconfiguration performance. For example, if the DGs do not affect the transient stability of the whole system significantly we can include an objective function similar to the switching cost (see Section-3.5) that incorporates the operating cost of DGs in reconfiguration problem. However, if large number of DGs are connected in the system then we need to consider their effect on transient stability of the whole system [33,19]. Similar to [19], we can consider the critical cleaning time (CCT) index of DG units as a reconfiguration factor to handle the effect of transient stability. In the future work, we will extend the proposed mathematical model to incorporate the effect of DGs in reconfiguration.

Conflict of interest

All authors declare no conflicts of interest in this paper.