Distribution system efficiency improvement using network reconfiguration and capacitor allocation

Electrical Power and Energy Systems 64 (2015) 457–468

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Distribution system efficiency improvement using network reconfiguration and capacitor allocation Hamid Reza Esmaeilian, Roohollah Fadaeinedjad ⇑ Department of Electrical and Computer Engineering, Graduate University of Advanced Technology, Kerman, Iran

a r t i c l e

i n f o

Article history: Received 5 August 2013 Received in revised form 15 June 2014 Accepted 16 June 2014

Keywords: Distribution network reconfiguration Capacitor allocation System efficiency Binary gravitational search algorithm

a b s t r a c t Distribution system utilities attempt to minimize cost with consideration of the system performance improvement. This paper employs the network reconfiguration and capacitor placement simultaneously to reduce energy losses and improve the system reliability subjected to satisfy operational and power quality constraints using a fuzzy approach. With respect to discrete and non-linear nature of the optimization problem, a Binary Gravitational Search Algorithm (BGSA) is utilized to solve the fuzzy multi-objective problem efficiently. The fast harmonic analysis method is used to perform harmonic power flow in the presence of shunt capacitors and non-linear loads. The state enumeration method based on Weibull–Markov stochastic model of the system components is adopted to assess reliability of different system configurations. Moreover, a new encoding strategy is proposed to boost the performance of network reconfiguration procedure. The IEEE 33-bus and an 83-bus practical distribution network of Taiwan Power Company (TPC) with a number of harmonic generating loads are utilized to test and validate the proposed method, and results are compared with Binary Particle Swarm Optimization (BPSO) algorithm and Binary Genetic Algorithm (BGA). Ó 2014 Elsevier Ltd. All rights reserved.

Introduction Deregulation pressure on Distribution Companies (DisCo) has led to cut cost and enhance the system efficiency in delivering power to customers. Among three power system sections (generation, transmission, and distribution), the most portion of the power loss (about 5–13% of the total power generation [1]) and the highest number of customer interruptions (approximately 80% of the total interruptions [2]) are related to distribution networks. Moreover, due to penetration of non-linear loads and voltage drop at end buses of radial feeders, power quality issues has been taken into account in the distribution network performance evaluation. Hence many alternatives such as network reconfiguration and capacitor placement are proposed to solve these problems. Distribution systems have been operated radially to facilitate their protection scheme and reduce the short circuit current. Therefore, each load point is fed by a route through the system components to the substation. So these systems have low reliability and high power loss. There are two types of switches in distribution networks, sectionalizing switches that are used to isolate the faulty network and tie switches are utilized to restore power ⇑ Corresponding author. Tel.: +98 342 6226517; fax: +98 342 6228018. E-mail addresses: [email protected] (H.R. Esmaeilian), fadaeind@yahoo. com (R. Fadaeinedjad). http://dx.doi.org/10.1016/j.ijepes.2014.06.051 0142-0615/Ó 2014 Elsevier Ltd. All rights reserved.

to customers through another route in order to reduce interruption duration in abnormal operation conditions. Since by status change of switches, the power flow to loads will be changed and consequently affects the power loss, voltages, harmonic distortion level, as well as the system reliability, hence in normal operation condition can improve the distribution network performance and reduce the cost by selecting the correct status of switches. The network reconfiguration was introduced by Merlin and Back [3] to reduce the power loss using branch-and-bound type heuristic technique. The switch exchange method was proposed by Civanlar et al. [4] and was extended by Goswami and Basu [5]. The algorithm was subsequently developed by Ababei and Kavasseri [6] to find sets of multiple switch exchanges that are implemented concurrently. These algorithms start with a radial configuration of the network and using heuristic formulas modify the network configuration in order to reduce the loss. Shirmohammadi and Hong [7] proposed a heuristic type algorithm to find the tie switch position in each loop to reduce the loss. The heuristic methods are usually fast but may not achieve the optimal configuration. Therefore, meta-heuristic algorithms have gradually been utilized to minimize the loss, such as GA [8,9], Ant Colony Algorithm (ACA) [1], Bacterial Foraging Optimization Algorithm (BFOA) [10], and Artificial Immune System (AIS) [11]. In different research, the network reconfiguration is applied as a multi-objective problem. In [12], a fuzzy multi-objective approach is considered for the power loss

458

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

reduction, node voltage deviation constraint, branch current loading constraint, and feeder load balancing. In [13,14], a network reconfiguration methodology is presented to maximize reliability of the power at the load points and minimize the system power loss. They considered models of various components to evaluate probabilistic reliability. In [15,16], a pareto-optimal microgenetic algorithm is employed to do a trade-off between the reliability indices and the power loss to find the optimal configuration. The network reconfiguration is an effective method to improve the system efficiency without any financial investment. However, reconfiguring network may not be able to satisfy the aimed loss reduction and power quality constraints. Therefore, the network reconfiguration is employed with capacitor placement or Distributed Generators (DG) to achieve better performance. In [17], the reconfiguration and DGs are used to minimize the power loss and balance load factor using ACA. In [18,19], the reconfiguration and capacitor placement are employed simultaneously to reduce the power loss. They showed that taking into account these two subjects simultaneously is more effective than considering them separately. A sensitivity index-based methodology is presented in [20] to minimize the energy loss using network reconfiguration combined with capacitor allocation. The shunt capacitors are applied for reactive power compensation locally to reduce the power loss and sustain the voltage profile within an acceptable operating limit. The shunt capacitors can be considered as redundant lines in the system inasmuch as reducing the reactive current leads to enhance load-carrying capability of the lines, as discussed in [21]. Therefore, they can also improve the system reliability. In a real distribution network, harmonic distortion levels are high due to harmonic injection by non-linear loads. If the harmonic sources in the system are not considered in capacitor placement, it may lead to an increase in harmonic distortion levels due to resonance between shunt capacitors and the various inductive elements in the system. In [22], optimal placement of shunt capacitor was applied to reduce the total cost of power loss and capacitor installation as well as the harmonic distortion level, using fast harmonic power flow method presented in [23]. Since reconfiguration alters the impedance structure of the network, so it affects the network harmonic distortion level. Nonetheless, in previous papers such as [18–20], the reconfiguration and capacitor placement were performed without harmonics consideration. In this paper, the network reconfiguration and capacitor placement are simultaneously employed to enhance the system efficiency in a fuzzy multi-objective optimization problem. The fuzzy membership functions are comprised of two categories: (1) minimizing the overall annual cost of the losses, expected outage cost (ECOST), and that of shunt capacitors, (2) satisfying operational and power quality constraints such as bus voltage limit and the system harmonic distortion level. The optimization problem is a combinatorial and non-differentiable problem, and due to discrete and non-linear nature, the BGSA is utilized to solve it. This paper is organized as follows: ‘Harmonic backward/ forward power flow’ presents the harmonic Backward/Forward power flow to calculate the network power loss and operating conditions. ‘Reliability assessment in radial distribution system’ presents the network reliability stochastic assessment to calculate the system reliability indices. The fuzzy multi-objective function is formulated in ‘Problem formulation’. ‘Application of binary gravitational search algorithm to solve reconfiguration and capacitor allocation problem’ describes the BGSA as an effective method to solve the optimization problem. Moreover, the agents encoding method is expressed in this section. ‘Test cases’ shows the simulation results of applying the proposed algorithm to two test systems. Furthermore, in this section the results are compared with other search algorithms. Eventually, the conclusion is outlined in ‘Conclusion’.

Harmonic backward/forward power flow In this paper, the harmonic analysis based on backward/forward sweep technique is used to perform power flow in distribution system in the presence of non-linear loads. Since the backward/forward sweep method converges very fast and needs low computational memory; it has widely been used in distribution networks to perform power flow. However, this method cannot be directly used for harmonic analysis because the harmonic currents absorbed by the shunt capacitors in network are unknown. In [23], the backward/forward sweep was developed to solve harmonic power flow problem. The proposed method is faster and more effective than other conventional harmonic analysis methods such as the admittance matrix. Because the procedure of the admittance matrix inverse is not necessary in this method. The algorithm can be divided into two parts, consisting of power flow calculation in fundamental harmonic and then in other harmonic orders. Fundamental power flow The algorithm starts with calculating the load current injection at all buses using Eq. (1). In backward sweep by direct application of KCL, the current of each branch can be calculated by summing current injected in the receiving bus and the set of the branch currents connected to the receiving bus, according to Eq. (2). Thereupon in forward sweep by direct application of KVL, the receiving bus voltage of each branch can be calculated by the difference between the sending bus voltage and voltage drop on the branch using Eq. (3). After performing these steps at each iteration, the mismatch between the calculated load power and the specified load power at all buses is calculated. If absolute value of the mismatch at all buses is less than certain tolerance value, then the iterations will stop, otherwise the steps will be repeated until the convergence is achieved.

Iki ¼

Si 

ð1Þ



V k1 i

ðkÞ

J L ¼ IðkÞ n þ

m X

ðkÞ

Jb

ð2Þ

V kn ¼ V kn1  Z L J kL

ð3Þ

b¼1

In the foregoing equations, V k1 denotes the voltage at node i at i ðkÞ the iteration k  1; Si is the specified power injected at node i; J L ðkÞ is the current of branch L at the kth iteration, In is the injected current at node n; m is the total branches connected to branch L on the receiving side, Z L is the series impedance of the branch L; n and n  1 are the receiving bus and the sending bus of the branch L, respectively. Fast harmonic power flow The algorithm steps are explained as follows: System harmonic current vector calculation The hth harmonic order current vector can be divided into two categories. The first category comprises the currents contributed by nonlinear and linear loads, ½Ih. The second category comprises the currents absorbed by shunt capacitors, ½Is. Therefore, the system harmonic currents vector at iteration k for the hth harmonic order can be expressed as:

2



I

 ðhÞ;k

3 ðhÞ;k Ih 6 7 ¼ 4 ... 5 IsðhÞ;k

ð4Þ

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

The nonlinear loads can be modeled as harmonic current sources injected into the system according to harmonic spectrum of the used device. The injection current of linear loads can be expressed as follows: ðhÞ;k Ihi

¼

V i ðhÞ;k

ð5Þ

ðhÞ

Zi

ðhÞ;k

where V i denotes the voltage of bus i at iteration k of the hth ðhÞ harmonic order, and Z i is the hth harmonic equivalent impedance of the linear load at bus i. The harmonic currents vector absorbed by shunt capacitors, ½Is, is obtained as the following equation:

h ih i    ðhÞ;k ðhÞ;k  HAsh ¼ HLF ðhÞ;k IsðhÞ;k Ih

The entry ai;j of matrix ½HAsh is the summation of the branches impedance in common path between capacitor i and linear or nonlinear load j toward slack bus. So instead of solving linear equations with dimension n, this method only needs to solve the linear equations with dimension s. Therefore, the expressed method is efficient. Backward sweep The branch currents can be calculated by the following equation.

h

ðhÞ;k

Bi;j

i

h iT    IðhÞ;k ¼ AðhÞ;k i;j

ð7Þ

ðhÞ

where ½Ai;j  is the coefficient vector of harmonic currents for branch between bus i and j at the hth harmonic order. If the harmonic current of a bus flows through the branch, then the corresponding position in ½Ai;j  will be 1, otherwise 0 is placed. Forward sweep The voltage drop in each branch can be calculated by the following equation. ðhÞ;k ðhÞ ðhÞ;k DV i;j ¼ Z i;j  Bi;j

ð8Þ

where Z i;j is the branch impedance between two buses i and j at the hth harmonic order. With respect to the system harmonic current vector, the bus voltages can be calculated by the forward sweep as follows:



V

ðhÞ;k





ðhÞ;k

¼ HA

  ðhÞ;k   I

voltages will be updated using Eq. (11), otherwise the algorithm will stop for the hth harmonic order. (Note: V h;1 ¼ 0 for all of buses.) ðhÞ;k DV ¼j V ðhÞ;kþ1  Vi j6 e i ðhÞ;kþ1

½V i

ð9Þ

ðhÞ;k

½HA  is a ðn  1Þ  m matrix, where n is the number of buses, and m is the total number of capacitors, linear loads, and nonlinear loads. This matrix represents the relationship between the bus voltage vector, ½V ðhÞ;k , and the system harmonic current vector, ½IðhÞ;k . Check for convergence After calculating the bus voltages, the difference between voltage at iteration k and k þ 1 for all of the network buses is computed. If the difference is bigger than the predetermined tolerance,

ðhÞ;k

 ¼ ½V i

ð10Þ

 þ DV

ð11Þ

After performing the power flow for all harmonic components, rms and Total Harmonic Distortion (THD) values of the bus voltages and the network power loss are calculated by the following equations.

V rmsi ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X ðhÞ j V i j2

ð12Þ

h¼1

ð6Þ

½HAsh is a s  n matrix and ½HLF is a s  s matrix that can be computed as follows; where s is the number of the installed shunt capacitors and n is the number of the system buses. 8 Diagonal entry ai;i : is the summation of the impedance of > > > > > capacitoriand the branches impedance in path > > > < capacitor i toward slack bus; ½HLF  ¼ > > > Off-diagonal entry ai;j : is the summation of the branches > > > impedance in common path between two > > : capacitors i and j toward slack bus:

459

THDi ð%Þ ¼ PLoss

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P ðhÞ 2 h¼2 j V i j ð1Þ

jV j X hi ¼ P1Loss þ PLoss

 100

ð13Þ ð14Þ

h¼2

Reliability assessment in radial distribution system The state enumeration method based on the Weibull–Markov stochastic model of the network components is used to evaluate the system reliability in this paper. The N  1 deterministic criterion does not distinguish between improbable events and more frequent events, hence may lead to undue expansion plans [24]. Since sequential Monte Carlo simulation is prohibitively timeconsuming, the contingency enumeration method is adopted to evaluate reliability of the different network configurations. The homogenous Markov models are unrealistic in the case of the repair duration modeling, and as a result the reliability worth indices cannot be accurately calculated by these models [25]. With this respect, the Weibull–Markov stochastic model is utilized to model failure frequency and repair duration of the system components states. In this paper, DIgSILENT Program Language (DPL) is linked to MATLAB in order to evaluate the system reliability. In such a way that the optimization algorithm determines which switches should be opened in reconfiguration problem, and through a file, DIgSILENT reads data and executes the reliability assessment for the considered configuration. The reliability of a network is evaluated by interruptions frequency and duration indices at the load points, and interruption cost to customers by utility. The customer-oriented reliability indices indicate the system performance in welfare of customers, that are weighted equally for all types of customer, such as F (systems average interruption frequency index) and T (system average interruption unavailability index) that are defined as follows [15]:

PNc i¼1 ki  kVAi F¼ P Nc i¼1 kVAi PN c i¼1 U i  kVAi T¼ P Nc i¼1 kVAi

ð15Þ ð16Þ

where ki and U i are average interruption frequency and average interruption time at the ith load point, respectively and Nc is the number of load points in the system. From the utility perspective, the interruption cost at the load points should be minimized. This requires an accurate calculation of the interruption costs. In many previous papers such as [15,21,26], the reliability assessment followed based on analytical or homogenous Markov methods. Since restoring power to the total loads is not possible by closing normally open switches (the loads located in downstream from failure in components which do not place in any loop), the repair duration of a system state is

460

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

stochastic; because the removing duration of the system failure is not the same in real conditions. Let’s consider component s with N states. Repair time of each state is characterized by shape factor (b) and scale factor (k), obtained from mean state duration (l) and state duration variance (r) [27], in Probability Density Function (PDF) of Weibull distribution. PDF for the repair time of component s can be expressed as follows [27]:

QN

bc;s N X eððT=kc;s Þ Þ c¼1 c;s ðTÞ= c;s PN c;s ðTÞ c¼1 c¼1 1= c;s

PrðU < TÞ ¼

l

c

c

l

ð17Þ

Fuzzy multi-objective formulation

ð18Þ

The objectives described in the previous section for the network reconfiguration and capacitor placement problem are mapped into fuzzy environment by membership functions in a single fuzzy satisfaction objective function as follows:

where;

cc;s ðTÞ ¼

"    bc;s !# kc;s 1 1 T C C ; bc;s bc;s kc;s bc;s

The interruption cost in load point i can be calculated as follows:

LPEICi ¼

M X

f k  Ei ðck ðdÞÞ

ð19Þ

k¼1

where f k is the frequency of the kth system state, M is the number of failures in ith load point, and Ei ðck ðdÞÞ is the expected value of the interruption costs that is calculated as follows:

Eðck ðdÞÞ ¼

Ns X PrðT j1 < Di < T j Þ  CðT j Þ

ð20Þ

j¼1

In Eq. (20), CðTÞ is Customer Damage Function (CDF) that is obtained by multiplying Sector Customer Damage Function (SCDF) in the average load in point i. The SCDF is defined based on the type of the customer that relates cost $/kW to interruption duration. For instance, SCDFs for customer types residential, commercial, and industrial are given in Table 1 [26]. Based on Table 1, CDF is defined as a vector of interruption durationsfT 1 ; . . . ; T Ns g, where N s is the number of intervals. The average load at the ith load point is calculated as Pi ¼ Peak Loadi  Load Factor [13]. Finally, the interruption cost paid to customers by utility (ECOST) is calculated by summing LPEICs at the load points. The shunt capacitors can enhance the system reliability indices by failure rate reduction of the network branches [21]. Flowing current through overhead lines and cables causes the resistive loss, and as a result leads to increase temperature. The high temperature will leads to life expectancy decrease of cables isolation, sag, and reducing the ground clearance of lines. Then by reducing reactive compoafter nent of the branch current can reduce the failure rate. If Ibefore andIr r are considered the reactive current of the ith branch before and after compensation by capacitors, then the compensated failure rate of the ith branch can be calculated as follows [21].

kcomp ¼ 0:85kuncomp þ 0:15kuncomp  i i i

Iafter r

ð21Þ

Ibefore r

function involves the overall annual cost of energy losses, customer interruption cost, and that of capacitors to be installed. The operational and power quality constraints such as THD and voltage profile of the network buses are considered to improve the network performance. The fuzzy set theory based on the decision-making method [28] is adopted in order to satisfy the objectives and soft constraints (instead of hard limits), the best compromise among them, as well as increasing efficiency of computation.

Problem formulation In this paper, the considered objectives to enhance the system efficiency are divided in two categories; the system cost minimization and the network performance improvement. The system cost

f ¼ max w1 lCost þ w2

1 min

20 min

60 min

120 min

180 min

Industrial Commercial Residential

1.625 0.381 0.001

3.868 2.969 0.093

9.085 8.552 0.482

14.44 16.141 1.959

19.8 23.73 3.437

! ð22Þ

subjected to the radial structure of distribution network in which all loads must be served. This constraint will be complied by encoding method expressed in ‘Application of binary gravitational search algorithm to solve reconfiguration and capacitor allocation problem’. In Eq. (22) lCost ; lV , and lTHD represent the membership functions for cost, voltage profile, and voltage THD, respectively. The wi coefficients are weighting factors that satisfy the condition: P3 i¼1 wi ¼ 1. The weighting coefficients are chosen by Decision Maker (DM) to find the solution that best satisfies the criteria [29,30]. Solving Eq. (22) yields a single solution that represents the best compromise considering the chosen weights. The proposed membership functions, on the contrary shape as vary from zero to one, are stated in the following sections. Membership function for cost minimization For the cost function a exponential membership function is defined as depicted in Fig. 1(a) and expressed as the following equation, where Costmin is the minimum system cost obtained in single objective optimization.

(

lCost ¼

1 exp



Costmin Cost Costmin



Cost 6 Cost min otherwise

ð23Þ

where;

Cost ¼ K E Eloss þ ECOST þ

nc X K ci  Q ci

ð24Þ

i¼1

In Eq. (24), K E is the energy cost per kW h, and Eloss is the network energy losses in a year that is calculated as Eloss ¼ 8760 loss load factor  P loss . The loss load factor can be calculated as G ¼ 0:5F þ 0:5  F 2 , where F is load factor [13]. K c is the capacitor annual cost per kVAr. Membership function for minimization of voltage constraint violation For the voltage constraints, it is appropriate to use of a trapezoidal membership function as shown in Fig. 1(b) and described in Eq. (25). V max and V min are the acceptable upper and lower bounds of the bus voltage that are considered 1 and 0.95 p.u., respectively. V  is the minimum bus voltage of the network embryonic situation, and V þ is the maximum predefined value of the bus voltage.

Table 1 SCDF for types of customer ($/kW) [26]. User sector

n 1X l þ w3 lTHD n i¼1 V i

lV i

8 0 > > > V i V  > > > > < V min V  ¼ 1 > > V þ V i > > > V þ V max > > : 0

Vi < V V  6 V i < V min V min 6 V i < V max V

max

6 Vi < Vþ

Vi > Vþ

ð25Þ

461

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

Fig. 1. Membership functions for (a) cost function, (b) rms voltage, (c) voltage THD.

Membership function for minimization of total voltage harmonic distortion The membership function ðlTHD Þ shown in Fig. 1(c), based on Eq. (26), is considered for THD. THDmax is the maximum acceptable limit (5% [22]), and THDþ is the maximum initial value of THD in the base case.

lTHD ¼

8 > <1 > :

THD < THDmax

THDþ THD THDþ THDmax

THDmax 6 THD < THDþ

0

THD > THDþ

ð26Þ

Load model There are different types of load in a real distribution network, such as residential, commercial and industrial. Since the load model type affects power flow equations and consumers interruption cost, and as a result on the network reconfiguration and capacitor allocation problem, hence, it should be considered in the optimization problem. In this paper, the static voltage dependency to load powers is modeled as polynomial model according to the following equations. The respective exponents of Eqs. (27) and (28) are shown in Table 2 [31].

a  bP  cP ! V P V V þ BP þ CP V0 V0 V0  aQ  bQ  cQ ! V V V AQ þ BQ þ CQ V0 V0 V0

P ¼ P0 AP Q ¼ Q0



ð27Þ ð28Þ

where the coefficients A, B, and C are percentages of load types connected to each bus, so that A þ B þ C ¼ 1 . Application of binary gravitational search algorithm to solve reconfiguration and capacitor allocation problem Overview of BGSA The GSA is a new meta-heuristic population search algorithm based on swarm intelligence, which was introduced by Rashedi et al. [32]. It is derived from the gravity and motion laws among

Table 2 Values of exponential parameters for the load types [31].

P Q

Industrial (a)

Commercial (b)

Residential (c)

0.1 0.6

0.6 2.5

1.7 2.6

masses as natural physical phenomenons. In this algorithm, the searcher agents are considered as masses in the search space that interact with each other based on the Newtonian laws, and their performance is evaluated by their masses amount. Based on the gravity law, ‘‘each particle attracts another by a gravity force that is directly proportional to the product of their masses and inversely proportional to the square of distance between them’’. The resultant force acting upon each agent from the others, leads to a movement toward the agents with heavier masses. The heavier mass corresponds to better solution of the problem. This relation is illustrated in Fig. 2. In GSA, each agent is specified by its position and mass. The position of the agent corresponds to a solution of the problem, and the mass value is determined using the fitness function so that bigger mass corresponds to a better solution. The heavier masses act the stronger gravity force upon other masses, so the agents tend to move toward them in order to search the space. Furthermore, the heavier masses have slower movement to search the space more locally. The GSA has many advantages than PSO [32], such as less memory consumption, considering fitness values and position of agents in updating procedure in order to achieve better solution. The main steps of GSA are as follows: Step Step Step Step Step Step

(1) (2) (3) (4) (5) (6)

Initialize the parameters and the agents position. Evaluate the fitness for each agent. Update GðtÞ, bestðtÞ and worstðtÞ of the population. Calculate the mass value for each agent. Update velocity and position of each agent. If the termination criterion has been met then end, otherwise go to step 2.

The process performs in a system with discrete time for each iteration of optimization algorithm. The algorithm steps are briefly described as follows.

462

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

law, the acceleration of agent i at time t and in direction d, adi ðtÞ, is calculated as follows:

adi ðtÞ ¼

Fig. 2. Each mass accelerates toward the resultant force from the other masses [32].

The position of the ith agent (mass) in a system with n dimension (the number of the problem variables) is defined as the following equation.

X i ¼ ðx1i ; . . . ; xdi ; . . . ; xni Þ for i ¼ 1; . . . ; N

ð29Þ

where xdi represents the position of agent i in dimension d, and N is the population size. After calculating the fitness value for each agent at time (iteration) tðfit i ðtÞÞ, the worst and best values in the population will be calculated as follows (in minimization problem):

bestðtÞ ¼ min fiti ðtÞ

ð30Þ

worstðtÞ ¼ max fiti ðtÞ

ð31Þ

i2f1;...;Ng

i2f1;...;Ng

The gravitational constant, G , is reduced with time to control the search accuracy, similar to temperature in the SA algorithm. It can be calculated by Eq. (32), where T is the number of iterations.

GðtÞ ¼ 2 

t T

ð32Þ

The value of the ith mass, M i ðtÞ , is calculated using the map of the fitness as the following equations:

fit i ðtÞ  worstðtÞ bestðtÞ  worstðtÞ q ðtÞ M i ðtÞ ¼ PN i j¼1 qj ðtÞ

qi ðtÞ ¼

ð33Þ ð34Þ

The force acting from mass j upon mass i, at time t and in direction dimension d, is defined as follows:

F di;j ðtÞ ¼ GðtÞ

 M i ðtÞ  M j ðtÞ  d xj ðtÞ  xdi ðtÞ Rij ðtÞ þ e

ð35Þ

In Eq. (35), e is a small constant and Rij ðtÞ is the Hamming distance between two agents i and j. All of forces acted on agent i in dimension d are weighted randomly in order to stochastic characteristic of the algorithm. The total force in dimension d can be obtained as Eq. (36).

F di ðtÞ ¼

X

randj  F di;j ðtÞ

ð36Þ

j2Kbest;j–i

where randj is a random number at interval [0, 1] and Kbest is the set of agents with the best fitness values (the biggest masses). In order to proper exploration and exploitation of the algorithm, the Kbest is time-varying as at the beginning of the algorithm it involves all agents, then it is decreased linearly with time. Using the motion

F di ðtÞ M i ðtÞ

ð37Þ

The position of a agent in every dimension of binary space is taken as 0 or 1. Therefore, the movement corresponds with value change from 0 to1 or vice versa. Hence, the velocity of each agent in every dimension, updated using probabilistic Eq. (38), is taken into a function to bound within interval [0, 1] as Eq. (39) [33]. Based on the hypercube framework, for small velocities, the probability of movement is near zero and for large velocities, the probability of the position change (xdi ) will go high. With respect to this point, the updated value of the ith agent position is calculated as Eq. (40). Moreover, to achieve a good converge rate, v di should be limited, means jv di j < V max (V max ¼ 6 [33]).

v di ðt þ 1Þ ¼ randi  v di ðtÞ þ adi ðtÞ Sðv di ðt þ 1ÞÞ ¼ j tanh v di ðt þ 1Þj

ð38Þ ð39Þ



if rand < Sðv di ðt þ 1ÞÞ then xdi ðt þ 1Þ ¼ complement xdi ðtÞ else xdi ðt þ 1Þ ¼ xdi ðtÞ

ð40Þ

Proposed approach The flowchart of the proposed approach to maximize the system efficiency by the network reconfiguration and capacitor placement simultaneously using BGSA is shown in Fig. 3. Considering close status for all switches, the network loops should be opened to maintain the radial topology of the network by opening a group of switches. Moreover, all of the network buses except the substation buses are candidate for capacitor placement. Hence, each agent is divided into two parts. First part specifies the number of switches that should be opened, and second part is considered to determinate the size of capacitors at each bus. Encoding strategy for the network reconfiguration based on graph theory is described in this section. Since the status of switches is randomly determined by meta-heuristic algorithms, hence it may lead to topologies of network that are not connective or tree. Fig. 4 shows some situations that lead to a non-radial topology of a typical distribution network. In many previous papers, nonfeasible agents have been removed by filters at each interaction of algorithm in order to avoid the evaluation of them. In large scale distribution systems, due to existence of complex loops with many common branches, the majority of the generated agents are not feasible and have to be removed of the population. Hence, it leads to a slow convergence rate. In [15], a feasible initial population is defined using the fundamental loops vectors of the network. Then new population will be created using ‘‘accentuated crossover’’ and ‘‘directed mutation’’ to avoid of nonfeasible agents. However, radiality of all topologies is not entirely provided, moreover, it requires to define a feasible initial population. In this paper, in order to speed up the computational time, nonfeasible agents in the population will be modified to feasible agents by applications of graph theory. In each mesh of the network should be opened a switch. Hence, each dimension of the agent is corresponded to a mesh of the network as it shows binary code related to switch number that is going to be opened in its corresponding mesh. Using this method, the number of the topologies that are not radial in the population will be decreased, because the branches which do not constitute loop with other branches have to be removed from the algorithm selection list. However, since each mesh has many common branches with other meshes, radiality of the network may not be

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

463

Fig. 3. Flowchart of applying network reconfiguration and capacitor allocation to maximize the system efficiency using BGSA.

Fig. 4. A typical distribution network.

satisfied. A subroutine is called at the beginning of the algorithm to determine the constituent branches of each mesh and graph adjacency matrix (adj) from the network connectivity matrix. The graph adjacency matrix is a n  n matrix, where n is the number of nodes. This matrix shows relationship between vertexes of graph. If there is an edge between the vertexes i and j, entry aij is 1, otherwise 0 is placed. Matrix Lmb specifies the constituent

branches of each mesh, where m is the number of meshes and b is the number of branches. If branch k is located in mesh i, entry Lik is 1, otherwise 0 is placed. From this matrix, the mesh common branches matrix is defined, ‘mb . This matrix expresses the common branches between each mesh with other meshes, S ‘i ¼ j¼1:m fLi \ Lj g. To assess the structures of the agents the followj–i

ing method is proposed.

464

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

Table 3 Simulation results of the 33-bus network using BGSA. Case

Base Case Case Case

case 1 2 3

Case 4

ENS (MW h/yr)

Ploss (kW)

Cost ($/yr)

THDmax%

Vmin (p.u.)

T (h/yr)

F (f/yr)

Tie switches

4.785 4.573 4.573 3.915

181.471 135.355 135.355 113.560

52068.568 43429.054 43429.054 37223.628

6.862 5.21 5.21 10.103

0.919 0.9418 0.9418 0.956

2.587 2.459 2.459 2.085

3.456 3.201 3.201 2.708

s33 ; s34 ; s35 ; s36 ; s37 s9 ; s14 ; s28 ; s32 ; s33 s9 ; s14 ; s28 ; s32 ; s33 s10 ; s14 ; s28 ; s32 ; s33

4.007

121.563

39248.639

4.967

0.9501

2.135

2.806

s9 ; s14 ; s28 ; s31 ; s33

Step (1) Determining the matrixes adj0 , L, and ‘ by the subroutine at the beginning of the algorithm. Step (2) choosing switches by the optimization algorithm (one in each mesh), S ¼ fs1 ; . . . ; si ; . . . ; sm g; si 2 Li ; si – sj . Step (3) Updating the adjacency matrix for the determined structure, adj1 . Step (4) Removing the single buses: if the nth row of matrix adj1 is zero, the branches connected to this bus will be specified (entries one in the n throw of matrix ½adj0  adj1 ). One of them is randomly chosen, for instance branch i. Then si is randomly changed in set S Li  sj 2S f‘i \ ‘j g, and ½adj1  will be updated.

Items

Case 1

BGA

Best Worst Average

The mentioned method is a simple method to implement so that experimentally with one or two changes in switches status, the network will be radial. This method ensures radiality of the network, and since the sets L and ‘ are specified at the beginning of the algorithm and also the nonfeasible agents are restructured at each iteration, hence this method speeds up the algorithm convergence. To encode the second part of agents for sizing of capacitors, a look up table is adopted. The sizes of the employed capacitors are as multiples of 150 kVAr (see Table 7 of Appendix), so that the capacitor size in each bus changes from 0 up to the last considered multiple. The number of dimensions in the second part of agents corresponds to the number of the network buses minus the substation buses.

1 0.99

Voltage (p.u.)

0.98 0.97 0.96 0.95 0.94 0.93 0.92 0.91 5

10

15

20

25

Bus Number Fig. 5. The IEEE 33-bus network bus voltage profiles.

30

kVAr

Bus

kVAr

– – – 6 30 3 6 22

– – – 450 900 300 750 300

– – – 12 29 4 9 –

– – – 300 600 750 150 –

BPSO

13.68 43429.054 44982.196 43853.743

Best Worst Average

43429.054 44272.449 43608.241

11.27

CPU times (min)

Case 3

0.9597 0.9467 0.9544 15.35

Best Worst Average

CPU times (min) BGSA

Case 2

43596.364 45282.246 44125.437

CPU times (min)

73.97

0.9600 0.9483 0.9564 13.02

9.43

Case 4

37292.387 37804.945 37501.134

37264.991 37654.133 37418.583 65.66

0.9600 0.9512 0.9583 10.88

37223.628 37517.498 37321.427 58.62

0.9719 0.9616 0.9694 85.45 0.9723 0.9643 0.9706 72.93 0.9725 0.9656 0.9714 64.13

1 0.9

Fitness Function

sj 2‘j

Bus

Table 4 Comparison of BGA, BPSO, and BGSA results of the 33-bus network.

sj 2‘j

Step (5) By technique expressed in [34], radiality of graph is evaluated. If the network is radial, in consequence the agent is feasible, otherwise go to next step and determine all isolated buses [34]. Step (6) Randomly choose one of the isolated buses that its degree in matrix ½adj0  adj1  is one. The connected branch to this bus in matrix ½adj0  adj1  is randomly S changed in set Li  sj 2S f‘i \ ‘j g, then go to step 5.

Capacitors

0.8 0.7 0.6 0.5 0.4 50

100

150

200

250

300

Iteration Fig. 6. Convergence characteristics of BGA, BPSO, and BGSA for case 4 of the 33-bus network.

In this paper, BGA and BPSO are applied to solve the problem besides BGSA. Similar to BGSA, these algorithms begin with an initial population of the problem candidate solutions and improve it towards the optimal solution in an iterative loop. In each iteration of BGA, the population members, called chromosomes, are selected based on the fitness evaluation and the new generation is produced by using the genetic operators, i.e., crossover and mutation. The BPSO versus BGA uses swarm intelligence to find the optimal solution. In each iteration of BPSO, the motion of a particle, as a member of the population, in the search space is affected by three factors including its previous motion, the best position visited by it thus far, and the best position visited in the entire population. Details of how these algorithms are implemented are described in [33]. The following representation shows how the problem decision variables, i.e. tie switches and size of capacitors, are encoded inside

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

a string in the algorithms, where si and cj denote the tie switch number of mesh i and the size of capacitor at bus j, respectively.

465

The polynomial load model is adopted in each bus according to ‘Load model’. The SCDF for types of customer is given in Table 1. The mean repair duration and the switching times are set to 2 h and 30 min, respectively. The load factor and energy price are considered 0.5 and 0.05 $/kW h, respectively. For the purpose of comparison, BGA and BPSO are applied to solve the problem too. The algorithms were implemented in MATLAB and DIgSILENT computing environment as the flowchart presented in Fig. 3, and run on a PC with 2.86 GHz and 1 GB RAM. Test case A

Test cases To demonstrate the effectiveness of the proposed method, it is tested on the IEEE 33-bus distribution network and an 83-bus practical radial distribution network of TPC in the presence of non-linear loads in several cases. The harmonic spectra of non-linear loads and annual cost of discrete steps of commercially available capacitor sizes are shown in Tables 7 and 8 of Appendix.

The IEEE 33-bus distribution network is studied as the first test case. The network data and the initial configuration of the network along with the protection and manoeuvre schemes are obtained from [15]. The total real and reactive power loads on the network are 3715 kW and 2300 kVAr, respectively. The coefficients in Eq. (27) are set to 30% residential, 40% commercial, and 30% industrial at each load point. The power factor ranges for types of customer are considered 0.95 6 for residential customers, 0.8 6 for

Fig. 7. The 83-bus distribution network of TPC.

466

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

Table 5 Simulation results of the 83-bus network using BGSA. Case

ENS (MW h/yr)

Ploss (kW)

Cost ($/yr)

Base case

8.081

491.4

118196.942

Case 1

7.556

441.391

Case 2

7.567

441.6

Case 3

6.592

321.077

83065.805

Case 4

6.677

321.363

83527.084

THDmax% (p.u.)

Vmin (h/yr)

T (f/yr)

F

Tie switches

Capacitors

8.01

0.934

0.579

0.853

107434.243

6.552

0.951

0.539

0.771

107522.711

5.287

0.951

0.539

0.772

s84 ; s85 ; s86 ; s87 ; s88 ; s89 ; s90 s91 ; s92 ; s93 ; s94 ; s95 ; s96 s7 ; s13 ; s34 ; s39 ; s42 ; s55 ; s62 s72 ; s86 ; s89 ; s90 ; s91 ; s92 s7 ; s13 ; s34 ; s38 ; s42 ; s55 ; s62 s72 ; s86 ; s89 ; s90 ; s91 ; s92

– – – – – –

– – – – – –

– – – – – –

– – – – – –

0.971

0.464

0.667

s7 ; s13 ; s34 ; s38 ; s55 ; s63 ; s72 s83 ; s86 ; s89 ; s90 ; s92 ; s95

13 31 79 12 14 41 71 19

1350 1800 1200 900 450 450 900 1350

28 55 7 22 64 18 69 –

1200 600 1500 900 1500 1050 1200 –

0.973

0.469

0.677

s7 ; s34 ; s38 ; s42 ; s55 ; s63 ; s72 s86 ; s88 ; s89 ; s90 ; s91 ; s92

1 7 71 12 14 32 34 81

1800 1500 1650 1200 600 450 150 1200

79 28 54 31 39 2 20 –

1200 900 1800 1200 600 300 1800 –

Bus

10.01

4.866

Table 6 Comparison of BGA, BPSO, and BGSA results of the 83-bus network. Items BGA

Best Worst Average

CPU times (min) BPSO

Best Worst Average

CPU times (min) BGSA CPU times (min)

Best Worst Average

kVAr

Bus

kVAr

Table 8 Harmonic spectra of non-Linear loads [35,36].

Case 1

Case 2

Case 3

Case 4

107598.436 110033.901 108232.393 53.32 107434.243 109725.385 107933.456 48.49 107434.243 109013.680 107689.932 44.81

0.9874 0.9631 0.9806 58.98 0.9875 0.9693 0.9825 53.14 0.9879 0.9773 0.9850 47.59

83123.141 83713.578 83369.726 261.54 83099.824 83606.626 83293.378 241.69 83065.805 83451.512 83185.430 225.63

0.9966 0.9829 0.9928 285.68 0.9969 0.9863 0.9940 259.92 0.9971 0.9892 0.9958 238.79

commercial customers, and 0.6 6 for industrial customers. To study the impact of harmonics on the network operation, the main harmonic producers in the system are considered as 15% fluorescent, 10% PC, and 10% Adjustable Speed Drive (ASD) in buses 14, 22, and 24; and 15% ASD in bus 30. The maximum capacitor size for sitting at each bus is set at 7  150 kVAr. Meanwhile, the total reactive power injection by capacitors should not exceed the total reactive power demand in the network. The studied cases are as follow: Case 1: the network reconfiguration to minimize cost. Case 2: the network reconfiguration to enhance system efficiency. Case 3: the network reconfiguration and capacitor placement to minimize cost (Eq. (24)). Case 4: the network reconfiguration and capacitor placement to enhance system efficiency (Eq. (22)). The weighting coefficients for cases 1 and 3 are considered as w1 ¼ 1 and w2 ¼ w3 ¼ 0; for case 2,

Harmonic order

non-linear loads ASD

1 5 7 11 13 17 19 23 25 29 31 35 37 41 43 47 49

Fluorescent

PC Workstation

Mag (%)

Phase (deg)

Mag (%)

Phase (deg)

Mag (%)

Phase (deg)

100 82.8 77.5 46.3 41.2 14.2 9.7 1.5 2.5 0 0 0 0 0 0 0 0

0 135 69 62 139 9 155 158 98 0 0 0 0 0 0 0 0

100 8.6 2.9 1.4 0.8 0.2 0.5 0.2 0 0 0 0.1 0 0 0 0.1 0

2 128 73 6 29 17 56 65 0 0 0 42 0 0 0 137 0

100 3.9 1.4 0.9 0.5 0.5 0.3 0.3 0.3 0.4 0.4 0.4 0.2 0.4 0.4 0.2 0.2

0 150 28 176 60 139 126 112 155 38 134 44 78 102 39 163 127

w1 ¼ 0:78; w2 ¼ 0:1, and w3 ¼ 0:12; and for case 4, w1 ¼ 0:52; w2 ¼ 0:18, and w3 ¼ 0:3. The computational results for the considered cases using BGSA are shown in Table 3. In case 1, the network power loss and ENS are reduced by 25.41% and 4.431% of the initial configuration, respectively. The network operational conditions including THD, voltage deviation, T, and F reliability indices are reduced to 24.07%, 28.15%, 4.948%, and 7.378% respectively. In this case, the obtained cost saving value is 8639.514 $/yr and the minimum bus voltage of

Table 7 Annual costs of shunt capacitors [22]. Q c (kVAr) kc ($/kVAr) Q c (kVAr) kc ($/kVAr)

150 0.5 1050 0.228

300 0.35 1200 0.17

450 0.253 1350 0.207

600 0.22 1500 0.201

750 0.276 1650 0.193

900 0.183 1800 0.187

467

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

the system is 0.9418 p.u., which occurs at bus 32. Meanwhile, the best result acquired in case 2 is the same as case 1. Therefore, an effective reconfiguration leads to the reduction of the cost and the improvement of the system performance. In order to achieve more appropriate operational conditions and more cost saving, the capacitor placement is employed in cooperation with the network reconfiguration. In case 3, THD is increased to 47.231% , although the loss and ENS are reduced. Therefore, this combination of switches and capacitors is not practical. In case 4 , the best results are obtained. The network power loss, ENS, THD, voltage deviation, T, and F are reduced to 33.01%, 16.26%, 27.62%, 38.4%, 17.47%, and 18.81%, respectively. In this case, the obtained cost saving is 12819.93 $/yr and the minimum bus voltage of the system is 0.9501 p.u., which occurs at bus 32. The network bus voltage profiles for cases: base, 2, and 4 are compared in Fig. 5. For the comparison, the obtained computational results of applying three algorithms BGA, BPSO, and BGSA for 30 runs are tabulated in Table 4, where the CPU time is the average running time when the best solution is achieved. The papulation size in each algorithm is set to 15 for the reconfiguration problem, and 30 for

the reconfiguration and capacitor placement. In BGA, one-point crossover rate 0.9, uniform mutation probability 0.05, and roulette wheel selection were adopted. In BPSO algorithm, c1 ¼ c2 ¼ 2 and the inertia factor is linearly decreased from 0.2 to 0.9. The convergence rates of the BGA, BPSO, and BGSA for case 4 of the 33-bus network are depicted in Fig. 6. It is clear that the BGSA obtains better results in shorter time compared to other methods. Test case B The second test case is an 11.4-kV 83-bus practical radial distribution network of TPC [1]. The initial configuration along with protection scheme is shown in Fig. 7. This network has 13 tie switches, 83 sectionalizing switches, and 11 feeders. The total active and reactive power loads on the network are 28,350 kW, 20,700 kVAr, respectively. The residential, commercial, and industrial portions at each load point are set to 30%, 40% and 30% for feeders 1–6, and 30%, 50% and 20% for feeders 7–11, respectively. The harmonic main producers in the network are considered as 15% ASD in buses 6, 75, and 81; 30% ASD, 10% PC, and 15% Fluorescent in buses 34,

Table 9 Line and load data of the 83-bus distribution network. Branch Sd.-Rv.

B1–1 1–2 2–3 3–4 4–5 5–6 6–7 7–8 7–9 7–10 B2–11 11–12 12–13 12–14 B3–15 15–16 16–17 17–18 18–19 19–20 20–21 21–22 21–23 23–24 B4–25 25–26 26–27 27–28 28–29 B5–30 30–31 31–32 32–33 33–34 34–35 35–36 36–37 37–38 38–39 39–40 38–41 41–42 B6–43 43–44 44–45 45–46 B7–47 47–48

R + jX (X)

0:1944 þ j0:6624 0:2096 þ j0:4304 0:2358 þ j0:4842 0:0917 þ j0:1883 0:2096 þ j0:4304 0:0393 þ j0:0807 0:0405 þ j0:138 0:1048 þ j0:2152 0:2358 þ j0:4842 0:1048 þ j0:2152 0:0786 þ j0:1614 0:3406 þ j0:6944 0:0262 þ j0:0538 0:0786 þ j0:1614 0:1134 þ j0:3864 0:0524 þ j0:1076 0:0524 þ j0:1076 0:1572 þ j0:3228 0:0393 þ j0:0807 0:1703 þ j0:3497 0:2358 þ j0:4842 0:1572 þ j0:3228 0:1965 þ j0:4035 0:131 þ j0:269 0:0567 þ j0:1932 0:1048 þ j0:2152 0:2489 þ j0:5111 0:0486 þ j0:1656 0:131 þ j0:269 0:1965 þ j0:396 0:131 þ j0:269 0:131 þ j0:269 0:0262 þ j0:0538 0:1703 þ j0:3497 0:0524 þ j0:1076 0:4978 þ j1:0222 0:0393 þ j0:0807 0:0393 þ j0:0807 0:0786 þ j0:1614 0:2096 þ j0:4304 0:1965 þ j0:4035 0:2096 þ j0:4304 0:0486 þ j0:1656 0:0393 þ j0:0807 0:131 þ j0:269 0:2358 þ j0:4842 0:243 þ j0:828 0:0655 þ j0:1345

k (f/yr)

0.1269 0.1368 0.1540 0.0599 0.1368 0.0257 0.0264 0.0685 0.1540 0.0685 0.0513 0.2224 0.0171 0.0513 0.0740 0.0342 0.0342 0.1026 0.0257 0.1112 0.1540 0.1026 0.1284 0.0855 0.0371 0.0685 0.1625 0.0318 0.0855 0.1284 0.0855 0.0855 0.0171 0.1112 0.0342 0.3251 0.0313 0.0313 0.0513 0.1368 0.1284 0.1368 0.0318 0.0257 0.0855 0.1540 0.1587 0.0428

Receiving node P (kW)

Q (kVAr)

0 100 300 350 220 1100 400 300 300 300 0 1200 800 700 0 300 500 700 1200 300 400 50 50 50 50 100 100 1800 200 0 1800 200 200 800 100 100 20 20 20 20 200 50 0 30 800 200 0 0

0 50 200 250 100 800 320 200 230 260 0 800 600 500 0 150 350 400 1000 300 350 20 20 10 30 60 70 1300 120 0 1600 150 100 600 60 60 10 10 10 10 160 30 0 20 700 150 0 0

Branch Sd.-Rv.

48–49 49–50 50–51 51–52 52–53 53–54 54–55 B8–56 56–57 57–58 58–59 59–60 60–61 61–62 62–63 63–64 B9–65 65–66 66–67 67–68 68–69 69–70 70–71 71–72 B10–73 73–74 74–75 75–76 B11–77 77–78 78–79 79–80 80–81 81–82 82–83 5–55 7–60 11–43 12–72 13–76 14–18 16–26 20–83 28–32 29–39 34–46 40–42 53–64

R + jX (X)

0:0655 þ j0:1345 0:0393 þ j0:0807 0:0786 þ j0:1614 0:0393 þ j0:0807 0:0786 þ j0:1614 0:0524 þ j0:1076 0:131 þ j0:269 0:2268 þ j0:7728 0:5371 þ j1:1029 0:0524 þ j0:1076 0:0405 þ j0:138 0:0393 þ j0:0807 0:0262 þ j0:0538 0:1048 þ j0:2152 0:2358 þ j0:4842 0:0243 þ j0:0828 0:0486 þ j0:1656 0:1703 þ j0:3497 0:1215 þ j0:414 0:2187 þ j0:7452 0:0486 þ j0:1656 0:0729 þ j0:2484 0:0567 þ j0:1932 0:0262 þ j0:0528 0:324 þ j1:104 0:0324 þ j0:1104 0:0567 þ j0:1932 0:0486 þ j0:1656 0:2511 þ j0:8556 0:1296 þ j0:4416 0:0486 þ j0:1656 0:131 þ j0:264 0:131 þ j0:264 0:0917 þ j0:1883 0:3144 þ j0:6456 0:131 þ j0:269 0:131 þ j0:269 0:131 þ j0:269 0:3406 þ j0:6994 0:4585 þ j0:9415 0:5371 þ j1:0824 0:0917 þ j0:1883 0:0786 þ j0:1614 0:0524 þ j0:1076 0:0786 þ j0:1614 0:0262 þ j0:0538 0:1965 þ j0:4035 0:0393 þ j0:0807

k (f/yr)

0.0428 0.0257 0.0513 0.0257 0.0513 0.0342 0.0855 0.1481 0.3508 0.0342 0.0264 0.0257 0.0171 0.0685 0.1540 0.0158 0.0318 0.1112 0.0793 0.1428 0.0318 0.0476 0.0371 0.0171 0.2115 0.0211 0.0371 0.0318 0.1640 0.0846 0.0318 0.0855 0.0855 0.0599 0.2053 0.0855 0.0855 0.0855 0.2224 0.2994 0.3508 0.0599 0.0513 0.0342 0.0513 0.0171 0.1284 0.0257

Receiving node P (kW)

Q (kVAr)

0 200 800 500 500 500 200 0 30 600 0 20 20 200 300 300 0 50 0 400 0 0 2000 200 0 0 1200 300 0 400 2000 200 500 100 400 – – – – – – – – – – – – –

0 160 600 300 350 300 80 0 20 420 0 10 10 130 240 200 0 30 0 360 0 0 1500 150 0 0 950 180 0 360 1300 140 360 30 360 – – – – – – – – – – – – –

468

H.R. Esmaeilian, R. Fadaeinedjad / Electrical Power and Energy Systems 64 (2015) 457–468

41, and 44. The maximum size of capacitor is considered 12  150 kVAr. The computational results of BGSA for the studied cases are given in Table 5. In case 1, the network power loss, ENS, THD, voltage deviation, T, and F indices are reduced to 10.177%, 6.497%, 18.2%, 25.757%, 6.91%, and 9.61%, respectively. In case 2, by replacing s38 with s39 to normally open, THD is reduced by 19.31% in comparison to case 1, although the cost is increased by negligible amount of 0.082%. Therefore, operation of the network is more desirable in this case. In case 3, the maximum cost saving is acquired, whereas THD is increased by 24.97% compared to the base case. Therefore, the operation of the network is not allowed in this case. The most satisfactory results are obtained in case 4 so that the network power loss, ENS, THD, voltage deviation, and T and F reliability indices are reduced to 34.602%, 17.374%, 39.25%, 59.09%, 18.99%, and 20.633%, respectively. In this case, the cost saving value is 34669.86 $/yr. To verify the performance of the proposed algorithm, the test case was repeatedly solved for 30 times through the BGA, BPSO, and BGSA algorithms. The computational results including the best, worst and average values of the objective function and the average CPU time for four cases of three algorithms are shown in Table 6. The population size for cases 1 and 2 is set to 20, and 40 for cases 3 and 4. The other parameters are the same as test case A. Comparing the obtained results by three algorithms, it is observable that BGSA obtains the better results than the others. Conclusion A fuzzy multi-objective model has been presented in this paper to reduce cost and improve efficiency of the distribution systems by applying the network reconfiguration and capacitor placement simultaneously. The cost function comprises the overall annual cost of the energy loss, ECOST, and shunt capacitors to be installed. This paper has shown that if the power quality constraints are neglected in the network operation, it may lead to undesirable operational conditions. The fast harmonic power flow and the state enumeration method based on Weibull–Markov stochastic model of the network components have been adopted to calculate the terms of the objective function. The BGSA as a new powerful swarm intelligence algorithm has been employed to solve the problem. The effectiveness of this algorithm has been demonstrated by comparing the obtained results of two test systems with BGA and BPSO. Moreover, an efficient method has been proposed to modify the nonfeasible agents’ structures to feasible agents in the network reconfiguration in order to shorten the computational time. Appendix Annual costs of shunt capacitors and harmonic spectrum of the non-linear loads employed in this paper, as well as the 83-bus distribution network data are shown in Tables 7–9. References [1] Abdelaziz A, Osama R, El-Khodary S. Reconfiguration of distribution systems for loss reduction using the hyper-cube ant colony optimisation algorithm. IET Gener Transm Distrib 2012;6(2):176–87. [2] Chowdhury A, Koval D. Power distribution system reliability. John Wiley & Sons; 2009. [3] Merlin A, Back H. Search for a minimal loss operating spanning tree configuration for an urban power distribution system. In: Proc of the power systems computation conf (PSCC); 1975. p. 1–18. [4] Civanlar S, Grainger J, Yin H, Lee S. Distribution feeder reconfiguration for loss reduction. IEEE Trans Power Deliver 1988;3(3):1217–23. [5] Goswami S, Basu S. A new algorithm for the reconfiguration of distribution feeders for loss minimization. IEEE Trans Power Deliver 1992;7(3):1482–91.

[6] Ababei C, Kavasseri R. Efficient network reconfiguration using minimum cost maximum flow-based branch exchanges and random walks-based loss estimations. IEEE Trans Power Syst 2011;26(1):30–7. [7] Shirmohammadi D, Hong H. Reconfiguration of electric distribution networks for resistive line loss reduction. IEEE Trans Power Deliver 1989;4(2):1492–8. [8] Mendoza J, Lopez R, Morales D, Lopez E, Dessante P, Moraga R. Minimal loss reconfiguration using genetic algorithms with restricted population and addressed operators: real applications. IEEE Trans Power Syst 2006;21(2):948–54. [9] Torres J, Guardado J, Rivas-Dvalos F, Maximov S, Melgoza E. A genetic algorithm based on the edge window decoder technique to optimize power distribution systems reconfiguration. Int J Electric Power Energy Syst 2013;45(1):28–34. [10] Kumar K, Jayabarathi T. Power system reconfiguration and loss minimization for an distribution systems using bacterial foraging optimization algorithm. Int J Electric Power Energy Syst 2012;36(1):13–7. [11] Oliveiraa L, Oliveiraa E, Gomesa F, Silva I, Marcatoa A, Resendeb P. Artificial immune systems applied to the reconfiguration of electrical power distribution networks for energy loss minimization. Int J Electric Power Energy Syst 2014;56(1):64–74. [12] Das D. Reconfiguration of distribution system using fuzzy multi-objective approach. Int J Electric Power Energy Syst 2006;28(5):331–8. [13] Amanulla B, Chakrabarti S, Singh S. Reconfiguration of power distribution systems considering reliability and power loss. IEEE Trans Power Deliver 2012;27(2):918–26. [14] Kavousi-Fard A, Akbari-Zadeh M. Reliability enhancement using optimal distribution feeder reconfiguration. Neurocomputing 2013;106(1):1–11. [15] Mendoza J, Lopez M, Coello C, Lopez E. Microgenetic multiobjective reconfiguration algorithm considering power losses and reliability indices for medium voltage distribution network. IET Gener Transm Distrib 2009;3(9):825–40. [16] Mazza A, Chicco G, Russo A. Optimal multi-objective distribution system reconfiguration with multi criteria decision making-based solution ranking and enhanced genetic operators. Int J Electric Power Energy Syst 2014;54(1):255–67. [17] Wu Y, Lee C, Liu L, Tsai S. Study of reconfiguration for the distribution system with distributed generators. IEEE Trans Power Deliver 2010;25(3):1678–85. [18] Farahani V, Vahidi B, Abyaneh H. Reconfiguration and capacitor placement simultaneously for energy loss reduction based on an improved reconfiguration method. IEEE Trans Power Syst 2011;2(5):621–31. [19] Chang C. Reconfiguration and capacitor placement for loss reduction of distribution systems by ant colony search algorithm. IEEE Trans Power Syst 2008;23(4):1747–55. [20] Oliveira L, Carneiro S, de Oliveira E, Pereira J, Silva I, Costai J. Optimal reconfiguration and capacitor allocation in radial distribution systems for energy losses minimization. Int J Electric Power Energy Syst 2010;32(8):840–8. [21] Etemadi A, Fotuhi-Firuzabad M. Distribution system reliability enhancement using optimal capacitor placement. IET Gener Transm Distrib 2008;2(5):621–31. [22] Eaja A, El-Hawary M. Optimal capacitor placement and sizing in unbalanced distribution systems with harmonics consideration using particle swarm optimization. IEEE Trans Power Deliver 2010;25(3):1734–41. [23] Teng J, Chang C. Backward/forward sweep-based harmonic analysis method for distribution systems. IEEE Trans Power Deliver 2007;22(3):1665–72. [24] Brown R. The electric power engineering handbook, power system planning (reliability). CRC Press; 2001. [25] Casteren J, Bollen M, Schmieg M. Reliability assessment in electrical power systems: the Weibull–Markov stochastic model. IEEE Trans Ind Appl 2000;36(3). [26] Bin Y, Xiu-li W, Zhao-hong B, Xi-fan W. Distribution network reconfiguration for reliability worth enhancement. In: Int conf on power syst technol; 2002. p. 2547–50. [27] Casteren J. Assessment of interruption costs in electric power systems using the Weibull–Markov model. Technical report 381. Chalmers Univ of Tech, Sweden; 2003. [28] Bellman R, Zadeh L. Decision making in a fuzzy environment. Manage Sci 1970;17(4):141–64. [29] Lee K, El-Sharkawi M. Modern heuristic optimization techniques: theory and applications to power systems. Wiley-IEEE press; 2008. [30] Berredo R, Ekel P, Martini J, Palhares R, Parreiras R, Pereira J. Decision making in fuzzy environment and multicriteria power engineering problems. Int J Electric Power Energy Syst 2011;33(3):623–32. [31] Dasi J. Power system analysis: short-circuit load flow and harmonics. CRC Press; 2011. [32] Rashedi E, Nezamabadi-pour H, Saryazdi S. Gsa: a gravitational search algorithm. Inform Sci 2009;13(179):2232–48. [33] Rashedi E, Nezamabadipour H, Saryazdi S. BGSA: binary gravitational search algorithm. Springer; 2009. [34] Bounova G; 2009. . [35] Ulinuha A, Masoum M, Islam S. Hybrid genetic-fuzzy algorithm for volt/var/ total harmonic distortion control of distribution systems with high penetration of non-linear loads. IET Gener Transm Distrib 2011;5(4):425–39. [36] Grady W. Understanding power system harmonics. Dept of Electrical & Computer Engineering, Univ of Texas at Austin; 2012.