1. Introduction[1]

1.1.    Background and literature review

In recent decades, distributed generations (DGs) have increasingly penetrated transmission and distribution systems. The use of wind turbines (WTs) in generating electricity could decrease greenhouse gas emissions and lead to technical benefits if their locations and sizes are optimally determined in terms of technical issues of the system. However, concerns, such as ensuring the continuous connection of WTs to the grid in case of power quality disturbances, are essential issues that should be considered [1, 2].

In the literature, there are new challenges introduced in different study fields, such as microgrids [3] as well as the distribution [4] and transmission systems [5], due to the inherently stochastic output power generated by WTs. However, there are some studies in the literature which have proposed approaches for responding to such challenges. For example, authors in [5] have presented a stochastic mathematical model for the optimal allocation of energy storage units in transmission systems to reduce wind power spillage and load curtailment while managing congestion and voltage deviation problems. Taking such challenges into account, the work presented in [6] introduces a new methodology for increasing wind power integration in deregulated electricity markets upon considering voltage stability assessments and wind power uncertainties. Another example of the works in this field is the one presented in [7], in which an intelligent coordinated volt/VAr optimization approach is proposed to maximize energy savings for active distribution systems penetrated with WTs.

Reconfiguration is often employed in distribution networks and microgrids [8] to improve operating conditions in terms of the reliability [9] and/or voltage stability of the system [10]. A systematic and thorough review of distribution network reconfiguration techniques is presented in [11] for mitigating unbalance of distribution networks. In this review, the authors review some of the more recent and significant methods of distribution network reconfiguration, DG placement, and their sizing aimed at minimizing power losses and improving the voltage profile.

The effect of DG penetration on reconfiguration in distribution systems is another issue studied by researchers. For example, the optimal multi-criterion reconfiguration of radial distribution systems with renewable solar and wind energy sources using the weight factor method has been presented while considering reliability in [12]. Minimizing power losses, improving the system's voltage profile and stability, as well as enhancing reliability are the main objectives of the problem defined in [12]. However, not only are the function normalization and selection of weight factors wrong, but the models of WTs and PVs have not been defined in the reliability and stability problems. In another study, Guo et al proposed a novel distribution system reconfiguration model that included DGs and ESSs to enhance the service reliability and benefits of distribution networks with DGs and ESSs [13]. The optimization problem in [13] is defined to minimize the sum of the customer interruption costs, operation costs of switches, and depreciation costs of DGs and ESSs. However, it is assumed in [14] that only one single fault could occur to one branch; thisis insufficient for reliability analysis of the system. Furthermore, the problem has been solved as a mixed-integer problem, and uncertainties in the outputs of distribution generations and energy storage systems have not been considered. Authors in [14] have attempted to find the optimal location and number of switches in a distribution system using an optimization problem with a variety of objective functions to improve reliability and reduce ENS and other system costs. However, the effects of renewable energy sources, such as WTs, have not been considered in [14]. A dynamic network reconfiguration has been presented in [15] for time-varying load profiles in the presence of DGs, in which objective functions are defined to minimize power losses, operation costs of the distributed generation, and the amount of ENS. In [16], by combining the benefits of rearranging and replacing components of the system without changing their reliability, an integrated method is proposed to improve the reliability of the reconfigurable system cost-effectively. In [17], a framework is presented for solving the problem of optimal distribution system reconfiguration aimed at improving reliability. Jose et al [18] presented a path-based modeling framework for the distribution feeder reconfiguration problem and demonstrated its ability to optimize network losses, reliability, or both. A modified multi-objective particle swarm optimization technique has been proposed in [19] for active distribution system planning by considering reconfiguration, renewable energy sources, and DSTATCOM. This technique aims to enhance voltage stability, reduce pollution, improve reliability, and maximize financial benefits.

The above review of the most significant studies done in the field of distribution system reconfiguration to reduce operational costs in the presence of DGs shows that this field is to some extent mature; however, there is still a need for adopting an approach that takes into account almost all practical aspects. To find a solution for the reliable and cost-effective operation of a real distribution system, examining the real impacts of the uncertainties of DGs on the system operation using innovative techniques is necessary. This is an issue which is addressed in this paper.

1.2. Motivations and contributions

For electrical power distribution system companies (DISCOs) that possess and manage distribution systems from an economic and technical point of view, issues, such as reducing electrical losses, improving feeder reliability, and reducing operational costs are of great importance. To this end, DISCOs exploit reconfiguration as the simplest and most cost-effective way to achieve these goals without adding extra equipment to the system. During the initial configuration of a distribution system, there are switching devices already installed at maneuver points, some of which are closed, and the rest are open. In cases where specific optimization purposes are needed, the statuses of these switches can change to form the configuration that corresponds to the optimal solution.

On the other hand, the ever-increasing penetration of WTs in distribution systems has made DISCOs face new challenges when assessing the reliability of their systems. This is because, in this case, there are uncertainties in load supplying paths and in the portion of the loads supplied by WTs. However, to the best of our knowledge, such challenges and their effects have not been analyzed in the literature until now. Among few studies that carried out initial assessments in this field, one can refer to [20] in which it is assumed that WTs supply their neighboring loads in a way similar to the islanded condition. According to [20], a distribution system should be considered as a system partitioned into islands. In this regard, due to uncertainties in the power generated by WTs, a distribution system may experience the lack of power needed for supplying loads; besides, the sum of the power generated by WTs and the power supplied by the upstream network can exceed the total demanded load. However, it is worth noting that these assumptions may greatly affect reliability indices and may not be consistent with what is happening.

To eliminate the above-mentioned shortcomings in the reliability assessment of distribution systems penetrated with WTs, a new methodology is presented in this paper, which uses minimal cut sets as a tool for calculating reliability assessment indices, such as ENS. In the new methodology, minimal cut sets are found based on the graph theory, using only the topology of the study system. Since each minimal cut set includes a set of equipment in the system, whose individual outages cause a specific load not to be supplied, the reliability assessment of the system is straightforward.

With the rapid integration of WTs into distribution systems, the main concern is how to effectively deal with the stochastic nature of their output power when making decisions on the operation of distribution systems. In this regard, one of the main challenges is the impact of WTs and their stochastic nature on reliability indices and energy losses. In brief, the main contributions of this paper are highlighted as follows:

• In this paper, a new concept named IRZ is introduced to overcome the challenge posed by the stochastic output of WTs, by splitting the network. IRZ is a useful concept that helps to calculate the minimal cut set by determining the power flow direction. To this end, a distribution system penetrated with WTs is partitioned into several IRZs, each of which being a radial feeder with only one power flow direction.
• A simple method is proposed for calculating the reliability index based on minimal cut sets, as obtained using spanning trees.
• An optimization model is designed aimed at optimally selecting the statuses of switches and WTs' locations, taking into account the minimization of ENS and energy losses as the objectives of the multi-objective optimization problem. This optimization model takes advantage of the way the optimization problem is formulated and the solver, i.e.  MOHSA, is used to find a solution.

The rest of the paper is organized as follows: The proposed problem formulation is introduced in Section 2; section 3 presents the MOHSA algorithm as the solver of the proposed optimization problem; simulations and numerical analysis are elaborated in Section 4, and conclusions are presented in Section 5.

2. Problem formulation

In distribution systems, reconfiguration is performed by closing and/or opening switches when the system is operating in its normal conditions. However, the method of selecting these switches depends on the aims of the desired optimization problem. In this section, a multi-objective optimization problem is formulated. This optimization problem is aimed to find the best configuration for the system along with the optimal size and location of WT are determined so that both the energy not supplied (ENS) and energy losses per year are minimized. In the following subsections, these objectives and the way they are formulated in a form suitable for the proposed approach are discussed.

2.1. Minimization of power losses

One of the main objectives of reconfiguring a distribution system is to reduce active power losses. Active power losses for a system with b lines can be calculated as follows [10]:

Where R_i and B_i are the resistance and the current flowing through the i^th line, respectively. According to (1), we need a direct load flow (DLF) analysis to find the branch current required for calculating active power losses. It is worth noting that several methods such as the backward-forward method, the impedance matrix calculation method, and DLF have been presented so far for load flow analysis in distribution systems. Since one of the aims of this paper is to find the best distribution system configuration in terms of energy loss minimization per year, it is necessary to use a load flow method that is both fast and accurate enough. Such characteristics can be achieved using DLF in which the relationships between line currents and injected currents are presented using a matrix called the bus current injections and the branch current (BIBC) matrix. To provide a clear picture of DLF performance, we consider a 7-bus system shown in Fig. 1. For this system, the relationships between bus current injections and branch currents can be defined as follows [10]:

Fig. 1. The line diagram of the 7-bus system

where I_i is the load current at bus i. Equation (2) can be rewritten as shown in (3) [10].

It is worth noting that the dimensions of B, BIBC, and I are m × 1, m × (n - 1), and (n - 1) × 1, respectively, where m and n are the number of the branches and buses of the system, respectively.

To form the matrix for interpreting the relationships between branch currents and bus voltages (BCBV), first, the graph theory is used to find the shortest path between a load bus in the system and the generating bus. In graph theory, a graph is made up of dots connected by lines where two dots can be connected by only one line [10]. Besides, a dot is known as a vertex, and a line is known as an edge. This indicates that a graph is the pair G = (V, E) of sets where the elements of V are the vertices, nodes, or points of graph G, and the elements of E are its edges or lines.

Considering the concept of the graph theory, a distribution system with a radial structure is a specific graph called a spanning tree. According to [10], "if G is a graph with n vertices, a spanning tree will be a connected sub-graph that uses all vertices of G with                n-1 edges".

In a graph, a cut set is a disconnecting set that breaks a path between a sink node and a source node; in contrast, a minimal cut set has no subsets. As mentioned in the introduction section, a minimal cut set is a set of equipment in the system, whose individual outages create a specific load to be left unsupplied. To find these cut sets, i.e. the shortest path between a load point and an energy supply point, a graph theory-based technique is proposed in this paper. It should be noted that one of the main features of the proposed methodology is its ability to systemically find paths between a load and a generating unit as well as the direction of the power flow. Besides, it is used to calculate the reliability index using only the configuration of the system. Another feature of this algorithm is its high speed in calculating the reliability index. This feature is very effective in reducing the time needed to  achieve the final solution because there are too many different configurations in the problem to be analyzed.

After finding the shortest path between each load bus and the generating bus, the voltage differences for the distribution system can be obtained as in (4) [10]:

Or, it can be rewritten as [10]:

where and BCBV are matrices with dimensions (n - 1) × 1 and (n - 1) × m, respectively. It should be noted that in (4), Z_ij is the impedance of the line connecting buses i and j. Equations (3) and (5) can be merged to achieve the relationships between the voltage of buses and currents injected into buses as in (6) [10].

where, dimensions of DLF are (n - 1) × (n - 1).

Equation (6) states that to calculate line currents, the magnitude and phase angle of voltages of all buses are needed, which are unknown. Therefore, to find the voltages and currents, a simple iteration-based procedure as the one defined in [10] is used in this paper.

2.2. Reliability assessment at load buses

In this paper, based on the graph theory described in Section 2.1 and according to the network topology, minimal cut sets in the distribution system are identified. These cut sets are then used for computing the reliability index. To do so, the same procedure used in (4) is employed. Considering the system shown in Fig. 1, the relationship between the bus connected to the upstream network and any desired load busesin terms of the reliability index can be obtained using (7).

In (7), Q_k represents the unreliability index of the load at the k^th bus, and Q_mn is the index representing the unreliability between buses m and n. In general, (7) can be rewritten as in (8).

where  and  are matrices with dimensions (n - 1) × 1 and (n - 1) × m, respectively.

The non-zero elements in (7) in each row of the matrix on the right side of this equation represent minimal cut sets of the network. For example, the second row shows that the minimal cut set for the load at bus 3 is the set consisting of lines B2 and B3. As it can be seen, matrices in (4) and (7) have the same dimension, with the same rows and columns for non-zero elements. The relationship between (4) and (7) could simplify the problem and can in turn speed up the solving process, which is effective in solving problems, such as the reconfiguration problem.

After calculating the unreliability index of buses using (8), the unreliability index of the whole system can be obtained via (9) [14].

where Q_DS is the reliability index of the whole distribution system, and N is the number of the loads in the system. It is worth noting that in a distribution system, a first-order cut set is a cut set by which only one path between the upstream network and any loads could be created. This cut set is very useful and can be used to calculate the unreliability between the upstream network and the i^th load (Q^I), which will be discussed in the following sections. In Fig. 2, the first-order failure rate of the i^th component in a cut set with N_I components is represented by C_i^I. Therefore, due to the outage of one or more components from the first-order cut set (Q^I), the unreliability between the upstream network and the load in Fig. 2 can be obtained using (10) [18].

Where P (C_i^I)is the probability of the outage of the i^th component and is the joint probability of the outage of i^th and j^th components.

Since probabilities of the outage of the components are mutually exclusive, (10) can be rewritten as follows:

Fig. 2. The first-order cut set between the upstream network and a load

2.2.1 The impact of a WT on the reliability index and power flows

Let us consider the 7-bus example system (shown in Fig. 1) and assume that a WT is installed at bus 6 shown in Fig. 3. The output active power of a WT is stochastic due to its dependence on the naturally stochastic wind speed. According to Fig. 3, it is assumed that the active power generated by the WT can supply the loads at buses 5 and 6, as well as part of the load at bus 4, as it is determined by coefficient µ. Therefore, given the assumptions made above, the sample distribution system is divided into two zones based on the power supplying paths shown in Fig. 4. Accordingly, each of these zones is called an independent reliability zone (IRZ).

Fig. 3. An example system with a WT installed at bus 6

There is another concept named 'deep point' that corresponds to the bus with its load supplied from two sources [8]. Since the power demanded in bus 4 is supplied from both the WT and the upstream network, this bus is a deep point, which indicates that the reliability of bus 4 depends on both supplying paths. As a result, the reliability of this deep point is restricted to the path with the lowest reliability index. This, in turn, has the advantage of increasing the reliability of bus 4 because the portion of the load, which is not supplied, decreases due to the existence of two supplying paths. Accordingly, as the system is transformed into two radial IRZs, it is possible to employ the load flow algorithm as well as the reliability assessment algorithm as described in sections 2.1 and 2.2.

Fig. 4. IRZs for the sample system

2.2.2 Reliability models of equipment in a distribution system

Probabilistic reliability models of different components in a minimal cut set are calculated based on the history of failure and maintenance. Accordingly, availability of these components can be obtained from (12) as follows [20]:

where A is the unit availability,  MTTF  is the mean failure time, and MTTR  is the mean repairing time of component i; also,  and are different types of outages and repair rates, including , , , and . Besides, and are failure and recovery rates, respectively;  and are maintenance and recovery rates, respectively; and  is the end-of-life rate for component i.

It is worth noting that  and are determined by the response of element outages, such as maintenance, active and passive failures, and the end-of-life rate of the element. Active and passive failure modes are two common failure modes in power systems, which determine reliability models for equipment in power systems.  However, according to [9], passive events do not cause the CBs to operate; thus, they do not affect other healthy components; nevertheless, active events cause protection breakers to operate, and several other healthy components are removed from service.

Fig. 5 shows reliability models of equipment other than wind turbines [9]. Besides, Figs. 5 (a) and   5 (b) show reliability models of bus bars as well as lines and switches, respectively. As can be seen, in the reliability models of lines and switches, failure modes are divided into active and passive modes. Finally, Fig. 5 (c) shows the common model of the generator, transformer, and CB, which includes the two failure modes of active and passive.

(a)

(b)

(c)

Fig. 5. Reliability models of equipment other than wind turbines: (a) buses, (b) lines and switches, and (c) generators, transformers, and circuit breakers

For the reliability modeling of a wind turbine, it should be first noted that according to [21], the output active power of a WT is related to the wind speed at the turbine installation site. This relationship can be defined as follows:

where  is the rated power of the WT,  is the wind speed at the site of the i^th turbine, and  and are cut-in and cut-out wind speeds. Moreover, wind speed is intrinsically stochastic and follows the Weibull distribution as in (14) and (15), according to [22].

On the other hand, a WT can operate in both up and down modes with the probabilities of p and 1 – p, respectively. Therefore, the total power generated by WTs at time t can be a function of the wind speed and the availability mode of each WT, as shown in (16) [22].

where  is the total generated active power,  is wind speed, and  is the vector representing the mode of each WT (if the WT is in the up mode, it will be 1, and if the WT is in the down mode, it will be 0 ). If the failure time for a WT is T_i, we  have:

Besides, the cumulative distribution function (CDF) for a wind farm with M wind turbines is defined as in (18) [22]:

In (18), obtains values as follows [22]:

where

It should be noted that in (20),  is the CDF of the wind speed. Therefore, considering (18) and (20), the CDF for the power generated by a WT can be obtained as follows [22]:

3. The proposed optimization problem

This paper aims to minimize unreliability at load points as well as losses in the system. However, for better visualization of the case studies, unreliability at load points is converted into ENS per year, and losses are converted into energy losses per year [23]. Accordingly, an index named the load loss factor (LLF) is needed, which is described in this section.

In the related literature, to determine and reduce losses in a distribution system, several general models have been introduced to develop a relationship between the LLF and the load factor (LF). In this paper, a relationship as presented in (22) is used [23].

Different studies show that when the hourly loading, the load factor, and the loss factor vary, variations in coefficient a are very small; thus, a can be assumed to be constant. At this stage, given the constant value of a, it is possible to obtain a model for losses of the system. In this paper, the non-linear relationship between the LF and the LLF as defined in (22) as well as the value of a (a = 0.5) have been adopted from [23].

In the following section, an example is provided to show the procedure of transforming reliability and power losses into ENS and energy losses, respectively. Let us consider the 7-bus system shown in Fig. 1 and assume that the load factor (LF) is 0.5. Besides, let the connected maximum load be 1000 kW. Hence, using (22), the LLF is calculated at 0.375. Moreover, the average load equals 1000 × LF = 500 kW. Therefore, since downtime refers to a time when a system fails to perform its primary function, the average downtime per year is obtained as Q_DS × 8760, in this paper. Therefore, ENS is calculated as follows [22]:

Energy losses, i.e. P_Ty, in a year is obtained using (24) [22] as follows:

The proposed approach aims to perform the reconfiguration and optimal installation of WTs so that both ENS and energy losses per year are minimized. Reconfiguration is done based on the spanning trees obtained using the graph theory as explained in Section 2.2. Maintaining the radial structure of the distribution system is one of the major constraints that should be considered when reconfiguring the system. Open switches should be selected so that the structure of the distribution system remains radial. In this paper, the Matroid method, based on the graph theory, is used to ensure that the system remains radial after reconfiguration [24]. Therefore, open switches for configuration of the distribution system are selected based on the Matroid method.

Thus, the proposed optimization problem includes two separate objectives that should be minimized simultaneously. To solve this problem, MOHSA is utilized in this paper for its ability to solve an optimization problem with separate objective functions as those defined in (23) and (24).

The harmony search algorithm is formulated based on a collection of musicians who play their musical instruments (population members) to achieve pleasing harmony (a globally optimal solution). In this context, the evaluation of any collections is performed based on the anesthetic standard. Accordingly, this algorithm has a simple concept, possesses few parameters, and is easy to implement, which turns it into a suitable choice for real-world problems, like traveling salesmen.

The details of MOHSA are explained in [25], in which it is demonstrated that MOHSA performs better than algorithms, such as the genetic algorithm and particle swarm optimization, in solving optimal reconfiguration problems. Figs. 6, 7, and 8 illustrate the flowchart of MOHSA integrated with the proposed method for solving the optimization problem as defined in this Section. Moreover, Fig. 8 shows the details of the blocks of the procedure for calculation of the unreliability index for each IRZ.

Finally, it should be noted that the input variables of the problem are locations of WTs and the system configuration. Hence, the global solution arrived at by the solver determines the optimal locations of WTs and the best configuration of the distribution system.

4. Simulation results

In this section, the efficacy of the proposed approach is demonstrated using two distribution systems. The first system is the IEEE 33-bus system which is a standard test system used by numerous papers in this field; besides, the second system is the Bijanabad distribution system which is used as a real system to show the way the proposed methodology can lead to better evaluations in real cases. It is worth noting that reaching the maximum iterations, in all simulations, is considered as the stopping criterion in MOHSA.

4.1. Case study 1: IEEE 33-bus distribution system

The proposed approach in this paper is applied to the IEEE 33-bus distribution system. Active and reactive powers in this system are 3715 kW and 2300 kVAR, respectively. Besides, active and reactive power losses are 202.7 kW and 139.2 kVAr, respectively. Based on [23], the load factor (LF) is assumed to be 0.5. Hence, using (22), the LSF is calculated at 0.375.

It is assumed that all lines of the system can be opened and closed by switches. It is  clear that distribution systems are operated in the radial mode; thus, only the first-order cut set is needed for calculating reliability at each point of the system. To compare and verify the results under normal operating conditions, the failure rate and maintenance of the equipment of the IEEE 33-bus distribution system given in [23] are employed.

As explained in Section 2.2.2, the active power generated by a WT can be calculated using (13). In this paper, the rated, cut-in, and cut-out speeds of a WT are assumed to be 15, 3, and 25 m/s, respectively. Besides, it is assumed that the rated power of a WT is 500 kW with the probability (p) being equal to 0.9 and the CDF as presented in (25).

Fig. 6. Flowchart of MOHSA algorithm

Fig. 7. Flowchart of the overall algorithm for calculating the unreliability of a distribution system

Fig. 8. The procedure of calculating unreliability index of the system

4.1.1 Scenario #1

In the first scenario, the problem is solved regardless of the effects of WTs. In other words, this scenario aims to find the best configuration by which the reliability of the system has been improved as much as possible. However, losses are minimized, in case no WT exists in the system.

In this scenario, the format of the solution vector in the HM matrix is a vector containing the decision variables and objective functions as defined in (27).

In (27), OSi is the ith open switch, and OFi is the jth objective function. In this case study, there are five open switches whose statuses  determine the configuration of the system and two objective functions, i.e. ENS and PTy per year, as explained in (23) and (24).

 The Pareto front composed of four solutions obtained using MOHSA in this scenario is illustrated in Fig. 9. Besides, the characteristics of these solutions and the basic solution are detailed in Table 1. As shown in Table 1, by comparing the basic solution with four other solutions, it can be observed that uncertainties in the system increase, while losses decrease by moving from S1 to S4. This indicates that the effect of the solution (configuration) on these two objectives is reversed concerning each other. Besides, the Pareto front suggests that S1 is the best solution among all solutions in terms of reliability, while S4 has been the best solution in terms of loss reduction. Therefore, upon moving from S1 to S4, ENS increases and losses decrease. Also, Table 1 shows that the proposed optimization algorithm has been effective in controlling ENS and energy losses per year in this scenario. This is realized by the appropriate selection of open and closed switches, which is done by the proposed algorithm. A comparison is made between the first column of Table 1, i.e. the basic solution with the basic solution achieved in [23], which is the best solution for the IEEE 33-bus system in the literature. As further proof, the comparison reveals that our reliability assessment methodology has led to a solution similar to the one in [23]. Therefore, this verifies the effectiveness of the proposed method in making unreliability assessments; in fact, this indicates that the proposed method can be used for reliability assessments in any other distribution systems with or without WTs.

Fig. 9. Pareto front obtained for the IEEE 33 bus system (Scenario 2)

Table (1): Solutions obtained by applying MOHSA to the IEEE 33-bus system (Scenario 1)

4.1.2 Scenario #2

In the second scenario, reconfiguration and optimal WT placement are performed to achieve the goals of the proposed optimization problem. In this scenario, only one WT with a capacity of 500 Kw is considered to be ptimally located. The format of the solution vector in the HM matrix is shown in (28). This solution vector is similar to the one given in (27) with the exception that another variable related to the location of the WT, i.e. , has been added to the solution.

Results of this scenario are provided in Fig. 10 and Table 2. As it can be seen, the installation of the WT in the system has resulted in reducing ENS and annual energy losses effectively. This is because the WT installed at bus 31 has changed the path of supplying electrical power and thus created a new IRZ. Therefore, due to the existence of these two IRZs, system reliability indices have improved, and energy losses in the system have decreased. Through enlightenment, the numerical comparison of ENS values and annual energy losses between the two scenarios is provided below:

Fig. 10. Pareto front obtained for the IEEE 33-bus system (Scenario 2)

Table (2): Solutions obtained by applying MOHSA to the IEEE 33-bus system (Scenario 2)