Fuzzy stochastic long-term model with consideration of uncertainties for deployment of distributed energy resources using interactive honey bee mating optimization

Abstract

This paper presents a novel modified interactive honey bee mating optimization (IHBMO) base fuzzy stochastic long-term approach for determining optimum location and size of distributed energy resources (DERs). The Monte Carlo simulation method is used to model the uncertainties associated with long-term load forecasting. A proper combination of several objectives is considered in the objective function. Reduction of loss and power purchased from the electricity market, loss reduction in peak load level and reduction in voltage deviation are considered simultaneously as the objective functions. First, these objectives are fuzzified and designed to be comparable with each other. Then, they are introduced into an IHBMO algorithm in order to obtain the solution which maximizes the value of integrated objective function. The output power of DERs is scheduled for each load level. An enhanced economic model is also proposed to justify investment on DER. An IEEE 30-bus radial distribution test system is used to illustrate the effectiveness of the proposed method.

Keywords

component / distributed energy resources / fuzzy optimization / loss reduction / interactive honey bee mating optimization (IHBMO) / voltage deviation reduction / stochastic programming

Introduction

The uncertainties associated with load forecasting and the unavailability of equipment affect the system operation and planning decisions. Applying a proper method for modeling these uncertainties in the planning phase can reduce the risk of the decisions as well as the stochastic cost of operation. Ignoring the uncertainties in planning process will lead to a high risk and render the stochastic saving gained by applying the non-optimal decisions.

In this paper, a novel methodology to solve the complicated problem of finding optimal location and size of distributed energy resources (DERs) is presented considering the uncertainties associated with load forecasting. In the proposed stochastic planning scheme, the stochastic characteristics of load growth are simulated using the Monte Carlo simulation method. Each possible system state is represented by a scenario. The scenario reduction technique is employed to decrease the number of created scenarios.

The honey bee mating optimization (HBMO) is a relatively new technique that has been empirically shown to perform well on many of these optimization problems [ 1– 3]. Unfortunately, the standard HBMO algorithm often converges to local optima, especially while handling problems with more local optima and heavier constraints. In other words, standard HBMO greatly depends on its parameter adjustments, and often suffers from the problem of being trapped in the local optima so as to be premature convergence. Therefore, some modification has been required for the standard HBMO algorithm to improve its performance.

To elude this deficiency, in this paper, a novel HBMO algorithm is proposed for solving the automatic generation control (AGC) problem. This paper presents the interactive strategy by considering the universal gravitation between the queen and drone bees for the standard HBMO algorithm to retrieve the disadvantages, which is called the interactive honey bee mating optimization (IHBMO). In other word, the IHBMO introduces the concept of universal gravitation into the consideration of the affection between drone bees and the queen bee within a honey bee colony. By assigning different values of the control parameter, the universal gravitation should be involved for the IHBMO when there are various quantities of drone bees and the single queen bee.

The restructuring of power systems has caused an increasing interest in DERs. Successful application and worldwide tendency to DERs have led to emergence of new technologies in this area. Moreover the increasing awareness of environmental issues has more motivated the application of DERs [ 4].

Many benefits are gained by placement of distributed energy resources (DERs), yet they may cause some troubles in operation of distribution systems if they are installed without thorough consideration. Therefore, special cares should be taken in locating and sizing of DERs. A wide range of benefits, from loss reduction to voltage profile improvement, can be gained by placement of DERs in distributed systems. Therefore, the realm of study of distribution systems is replete with the works on solving the problem of DER placement with different objective functions. In Ref. [ 5] the most important benefits of DER have been modeled in economic terms. A set of indices were proposed in Ref. [ 6] for modeling and quantifying of the technical benefits of DERs.

Previous work on DER placement

A distributed generation (DG) capacity investment planning algorithm has been proposed in Ref. [ 7] using a new heuristic approach from the perspective of a distribution company. An optimal solution for DG capacity placement has been obtained through a cost-benefit analysis approach based on this new optimization model. The model aims to minimize the disco’s operation and investment costs of the disco as well as the cost of power loss, but the other benefits gained by DG placement are not considered. The analytical optimization process for determining the optimal location and size of DER has been aimed at minimization of power loss of distribution systems in Ref. [ 8]. Both radial and meshed distribution systems have been considered.

Cost-benefit analysis is one of the other approaches used in the literature. For example, a new heuristic approach for DG capacity investment planning from the perspective of a distribution company has been proposed in Ref. [ 9]. To solve the placement and sizing problem, a multi-period AC optimal power flow (OPF) has been proposed in Ref. [ 10]. The minimization of power loss is again the aim of the optimization algorithm. Loss reduction and reliability improvement have been also handled in Ref. [ 9] as a cost/worth analysis.

Two multi-objective formulations based on the genetic algorithm (GA) and a ϵ-constrained method as optimization techniques have been proposed in Ref. [ 11] for the placement and sizing of DER in distribution networks. The optimization process was a compromise between reduction in power losses, reliability improvement, and reduction in power to be purchased from the power market and minimization of the cost of network upgrading.

Facility-location problems have been studied in Ref. [ 12] while taking into account a hybrid uncertain environment involving both randomness and fuzziness. Since the fuzzy parameters of the locating problem are represented in the form of continuous fuzzy variables, the determination of value-at-risk is inherently an infinite-dimensional optimization problem that is not possible to be solved analytically. Therefore, a two-stage fuzzy facility location problem with value-at-risk has been proposed in Ref. [ 13], which results in a two-stage fuzzy zero-one integer programming problem. Both the costs and demands have been skillfully assumed to be fuzzy random variables in Ref. [ 14], a value-at-risk based fuzzy random facility location model has been built and a hybrid modified IHBMO approach has been proposed to solve such a complicated problem.

A fuzzified multi-objective GA based algorithm has been proposed in Ref. [ 15] for capacitor placement. Though the objective is to find the best location of capacitors in distribution networks, the model of objective function and the methodology can be used in DG placement problem. A method for reliability improvement and loss reduction has been suggested in Ref. [ 16] by installing fixed capacitor in a distribution system. Though the problem is the capacitor placement, the reduction of power loss at peak load was skillfully modeled as one of the objective functions.

Two types of load uncertainties for planning studies can be identified in the planning power systems, uncertainties associated with load forecasting and short-term uncertainties related to time/weather factors. Both of them have been considered in Ref. [ 17] and a simple GA-based optimization algorithm has been used to extract the best location and size of DERs in a distribution system. The Monte-Carlo simulation has been used to model the stochastic nature of the system loads. The stochastic nature of system components and load growth forecasting has been simulated and each possible system state has been represented by a scenario in Ref. [ 18]. A scenario reduction technique is used to decrease the computational burden of large number of scenarios.

Motivations and contributions of the present work

(1) A long-term stochastic load model was developed for DER placement considering stochastic load growth. A long-term stochastic model for system uncertainties is presented in this paper that is suited for application along with the IHBMO algorithm. The results of case studies show the necessity of stochastic modeling of the problem. Some other studies in the literature have considered the stochastic nature of the load and system components, but uncertainties are modeled in just one hour or just one year. In planning problems, it is necessary to model the uncertainties in the entire planning horizon.

(2) System operation was considered in planning phase for different system states. In this paper, the output power is scheduled for each load level to avoid the inconvenient rejection of more optimal solutions. In contrast, previous works considered the output power of DERs to be fixed at the maximum rated value while the load varied at each bus. This may render some optimal solutions infeasible due to violation of some constraints such as voltage magnitude limits in some load levels while in most of the other load levels there is no violation.

(3) Fuzzy optimization approach was applied to satisfy different objectives simultaneously in DER placement. So many studies have been conducted to reduce the cost of loss in distribution systems. The reduction of voltage deviation in order to reach a more flat voltage profile is also the subject of many studies in distribution systems. In this paper, the reduction of loss and power purchased from the electricity market, loss reduction in peak load level and reduction in voltage deviation are considered simultaneously as the objective functions. These objectives are first fuzzified and then integrated and introduced into an IHBMO algorithm in order to obtain the solution which minimizes the value of integrated objective function.

(4) Adoptive membership functions. The fuzzy approach has been applied in previous works such as Ref. [ 16] (for capacitor placement), but membership functions have been predefined. In this paper, a method is presented to find the appropriate membership functions in fuzzification process of objective functions. A method is also presented in order to make these objectives comparable with each other.

(5) Economic modeling was improved. Profit maximization is considered as one of the objective functions in order to justify the investment on DER installation comparing to the other investment opportunities, the benefit to cost ratio (BCR) is considered as a constraint whose value should be greater than a predefined value. This predefined value should be calculated based on the other investment opportunities.

The proposed method is tested on IEEE-30 bus radial distribution test system. The simulation results show the effectiveness of the proposed method in DER planning problems and the necessity of stochastic modeling.

Honey bee mating optimization

Standard HBMO

The concept of the original HBMO was provided in Ref. [ 18]. At the start of the flight, the queen is initialized with some energy content and returns to her nest when her energy is within some threshold from zero or when her spermatheca is full. In developing the algorithm, the functionality of workers is restricted to brood care, and therefore, each worker may be represented as a heuristic which acts to improve and/or take care of a set of broods. A drone mates with a queen probabilistically using an annealing function as

(1) \mathrm{Prob}\left(Q,D\right)={\mathrm{e}}^{-\frac{\mathrm{\Delta }\left(f\right)}{S\left(t\right)}},

where Prob (Q, D) is the probability of adding the sperm of drone D to the spermatheca of queen Q (that is, the probability of a successful mating); ∆( f ) is the absolute difference between the fitness of D (i.e., f (D)) and the fitness of Q (i.e., f (Q)); and S(t) is the speed of the queen at time t. It is apparent that this function acts as an annealing function, where the probability of mating is high when both the queen is still in the start of her mating-flight and therefore her speed is high, or when the fitness of the drone is as good as the queen’s. After each transition in space, the queen’s speed, S(t), and energy, E(t), decay using

(2) S\left(t+\mathrm{1}\right)={\alpha }_{\mathrm{H}\mathrm{B}\mathrm{M}\mathrm{O}}\times S\left(t\right),

(3) E\left(t+\mathrm{1}\right)=E\left(t\right)-{\mathit{\gamma }}_{\mathrm{H}\mathrm{B}\mathrm{M}\mathrm{O}},

where α_HBMO is the speed reduction factor and γ_HBMO is the amount of energy reduced after each transition (α, γ ∈ [0,1]).

Thus, the HBMO algorithm may be constructed with the following five main stages:

Step 1 The algorithm starts with the mating flight, where a queen (best solution) selects drones probabilistically to form the Spermatheca (list of drones). A drone is then selected from the list at random for the creation of broods.

Step 2 Creation of new broods by crossoverring the drones’ genotypes with the queen’s (breeding process). The breeding process can transfer the genes of drones and the queen to the jth individual based on Eq. (4).

(4) \mathrm{Child}={\mathrm{Parent}}_i+{\beta }_{\mathrm{HBMO}}\left({\mathrm{Parent}}_k-{\mathrm{Parent}}_i\right),

where β_HBMO is the decreasing factor, {\beta }_{\mathrm{HBMO}}\in \left(\mathrm{0},\mathrm{1}\right) .

Step 3 Use of workers (heuristics) to conduct local search on broods (trial solutions).

Step 4 Adaptation of workers’ fitness based on the amount of improvement achieved on broods as follows:

(5) {\mathrm{Brood}}_i^k={\mathrm{Brood}}_i^k\pm \left({\delta }_{\mathrm{HBMO}}+{ϵ}_{\mathrm{HBMO}}\right){\mathrm{Brood}}_i^k, {\delta }_{\mathrm{HBMO}}\in \left[\mathrm{0},\mathrm{1}\right],{ϵ}_{\mathrm{HBMO}}\in \left(\mathrm{0},\mathrm{1}\right).

Step 5 Replacement of weaker queens by fitter broods.

Interactive honey bee mating optimization

The standard HBMO technique has a flexible and well-balanced mechanism to enhance the global and local exploration abilities. The main disadvantage of the HBMO algorithm is the fact that it may miss the optimum, not strong enough to maximize the exploitation capacity and provide a near optimum solution in a limited runtime period. Nonetheless, the standard HBMO method is prosperous for finding the best answer in optimization problem. However, the algorithm makes use of an independent random, β_HBMO~U (0, 1). This factor used to affect the stochastic nature of the standard HBMO algorithm, as shown in Eq. (9). In other words, the original breeding process design of the drone bee’s movement only considers the relation between the queen bee, which is selected by the mating wheel selection, and the one selected randomly. Hence, it is not strong enough to maximize the exploitation capacity. To overcome this drawback, the IHBMO algorithm is proposed based on the structure of the original HBMO algorithm. In the other word, by employing the Newtonian law of universal gravitation [ 17] described in Eq. (11), the universal gravitations between the queen bee and the selected drone bees are exploited. The gravitational force between two particles (F₁₂) with the mass of the first and second particles (m₁ and m₂), and the distance between them (r₁₂) can be expressed as

(6) F_{\mathrm{12}}=G{\displaystyle \frac{m_{\mathrm{1}}{m}_{\mathrm{2}}}{r_{\mathrm{21}}^{\mathrm{2}}}}{\widehat{r}}_{\mathrm{21}},

where G is gravitational constant. {\widehat{r}}_{\mathrm{21}} denotes the unit vector in Eq. (7).

(7) {\widehat{r}}_{\mathrm{21}}={\displaystyle \frac{r_{\mathrm{2}}-r_{\mathrm{1}}}{|r_{\mathrm{2}}-r_{\mathrm{1}}|}}.

In the optimization procedure with the IHBMO algorithm, the mass m₁ is replaced by the symbol, F (Parent_i), which is the fitness value of the queen bee that picked by applying the mating wheel selection. The mass, m₂, is replaced by the fitness value of the randomly selected drone bee and is denoted by the parameter, F (Parent_k). The universal gravitation in Eq. (16) is formed in the vector format. Therefore, the quantities of it on different dimensions can be considered separately. Hence, r₂₁ is calculated by taking the difference between the objects only on the concerned dimension currently while the universal gravitation on each dimension is calculated separately. In other words, the intensity of the gravitation on different dimensions is calculated one by one. Therefore, the gravitation on the jth dimensions between parent_i and parent_k can be expressed in Eq. (8). Finally, Eq. (4) can be modified in iteration interval t with Eq. (9).

(8) F_{i{k}_j}=G{\displaystyle \frac{F\left({\mathrm{Parent}}_i\right)\times F\left({\mathrm{Parent}}_k\right)}{{\left({\mathrm{Parent}}_{kj}-{\mathrm{Parent}}_{ij}\right)}^{\mathrm{2}}}}·{\displaystyle \frac{{\mathrm{Parent}}_{kj}-{\mathrm{Parent}}_{ij}}{\left|{\mathrm{Parent}}_{kj}-{\mathrm{Parent}}_{ij}\right|}},

(9) {\mathrm{Child}}_{ij}\left(t+\mathrm{1}\right)={\mathrm{Parent}}_{ij}\left(t\right)+F_{i{k}_j}\left[{\mathrm{Parent}}_{ij}\left(t\right)-{\mathrm{Parent}}_{kj}\left(t\right)\right].

Extending the consideration of the universal gravitation between the drone bees, which is picked by the queen bee, and more than one drone bee is achievable by adding different F_ik [ Parent_k – Parent_i] into Eq. (9), the gravitation F_ik plays the role of a weight factor controlling the specific weight of [Parent_k – Parent_i]. The normalization process is taken in order to ensure that F_ik~U (0, 1). Through the normalization of F_ik, the gravitational constant (G) can be deleted.

In general, the whole process of the IHBMO algorithm can be summarized at the five main steps as follows:

(1) Generate an initial population. In this step, an initial population based on state variable is generated, randomly, the formulation of which is

(10) \mathit{D}=\left[{\mathit{D}}_{\mathrm{1}},{\mathit{D}}_{\mathrm{2}},\mathrm{...},{\mathit{D}}_n\right]{\mathit{D}}_i=\left(d_i^{\mathrm{1}},d_i^{\mathrm{2}},\mathrm{...},d_i^m\right),

where d_i^j is the jth state variable value of ith drone bee. For each individual (D_i) the objective function values are evaluated. The queen is chosen according to the best solution (minimum objective function value). The other solutions which are generated during this phase became the drones to be used during the mating flight (trial solutions). The speed of queen at the start and the end of a mating flight is generated in this step.

(2) Flight mating. At the start of the flight mating, the queen flies with her maximum speed and energy. A drone is selected from the population of drones with the largest Prob (Q,D). The mating probability is calculated based on the objective function values of the queen and the selected drone. If the mating is successful (i.e., the drone passes the probabilistic decision rule), the drone’s sperm is stored in the queen’s spermatheca and the queen’s speed and energy is reduced by Eqs. (2) and (3). Otherwise, the flight mating is continued until the speed S(t) is less than a threshold d or the queen’s spermatheca gets full.

(3) Breeding process. In this step, a population of broods is generated based on the mating between the queen bee and the drone bees stored in the queen’s spermatheca. The generation of new broods occurs by crossoverring the drones’ genotypes with the queen’s. The jth brood is generated by using Eq. (9).

(4) Adaptation of worker’s fitness. In this step, after new broods are generated, an attempt is made to improve both the new solutions and the best solution by using of a mutation procedure as Eq. (10).

(5) Replacement of weaker queens by fitter broods and check the termination criteria. If the new brood is better than the current queen, it takes the place of the queen. If the new brood fails to replace the queen, then this brood will be one of the drones in the next flight mating of the queen. Meanwhile, if the termination criterion is satisfied, then finish the algorithm. Otherwise, it discards all of the previous trial solutions (brood set). Then the generation of new drones is set and goes back to Step 2.

Stochastic long-term model

In the proposed stochastic planning model, each possible system state is called a scenario. These scenarios are created by the Monte Carlo simulation method to model long-term stochastic characteristics of the system components and bus loadings. The scenario reduction technique is used to decrease the number of created scenarios.

Monte Carlo simulation method

Figure 1 depicts the annual load duration curve (LDC) that is modeled as multiple load blocks. This model is used as an infrastructure in consideration of forced outages of lines and load forecasting inaccuracies.

Hours with similar loads are shown in each load block in Fig. 1. Future annual peak load and energy demand growths are equal to the base year values times the regarding growth rates [ 15]. In this study the growth rate is expressed as an average growth rate, denoted by AGRP and AGRE, and a random component, RCP_t and RCE_t for annual peak loads and total energy demands, respectively. Normally distributed random components with certain standard deviations are aggregated with the average growth rates to reflect the uncertainties in economic growth and/or weather changes. Random trajectories in the sth scenario and the yrth year, represented by P_yr,s for peak load, and represented by E_yr,s for total energy demand, are expressed in the Monte Carlo simulation based on [ 18] as

(11) P_{yr,s}=P_{\left(yr-\mathrm{1}\right),s}\left(\mathrm{1}+\mathrm{AGRP}+{\mathrm{RCP}}_{yr,s}\right),

(12) E_{yr,s}=E_{\left(yr-\mathrm{1}\right),s}\left(\mathrm{1}+\mathrm{AGRE}+{\mathrm{RCE}}_{yr,s}\right).

The future load block in scenario s and year yr, denoted by bl_yr,s, is calculated via a linear transformation of the base year load and is formulated as

(13) {\mathrm{bl}}_{yr,s}=a_s{\mathrm{bl}}_{\mathrm{0}}+b_s,

where

(14) a_s={\displaystyle \frac{E_{yr,s}-\mathrm{8760}{P}_{yr,s}}{E_{\mathrm{0}}-\mathrm{8760}{P}_{\mathrm{0}}}},

(15) b_s={\displaystyle \frac{P_{yr,s}{E}_{\mathrm{0}}-P_{\mathrm{0}}{E}_{yr,s}}{E_{\mathrm{0}}-\mathrm{8760}{P}_{\mathrm{0}}}}.

The bus load of bus z at each load block is calculated by multiplying the load distribution factor of that bus and load at each block,

(16) {\mathrm{PD}}_{z,b,yr,s}=D_{z,b,yr}{\mathrm{bl}}_{yr,s}.

Transmission line availability of line k at load block b in year yr denoted by UY_k,b,yr is used in the Monte Carlo simulation in which UY_k,b,yr= 1 indicates that the transmission line k is available at load block b in year yr while UY_k,b,yr= 0 indicates otherwise. Consequently, a scenario is consisted of RCP_yr,s, RCE_yr,s and UY_k,b,yr,s.

Scenario generation

For each uncertain variable, different states are considered, each with a corresponding occurrence probability. Figure 2 shows these states for system demand (as an example).

The sampling process in a non-sequential Monte Carlo is an iterative process, which is performed by generating a vector of random values between 0 and 1, each regarding to one of the uncertain variables. This vector is used to draw a sample. For example the value of demand (ith random variable) would be equal to \mu +st\cdot \sigma if Eq. (10) is fulfilled.

Modeling the load uncertainties using a normal distribution causes some errors and difficulties. In practical power systems, historical data on bus loads is available. The load historical samples at each bus in several groups can be categorized based on their differences as a measure. The mean value of each group can be considered as different states. Occurrence probability of each group can be defined as the number of historical samples in this group divided by the total number of samples. However, this method of modeling cannot be used in the case studies in this paper, as a result of the lack of historical data. Therefore, a normal distribution is applied to model continuous uncertainties.

There are different approaches for modeling the uncertainties associated with the un-availabilities of equipments. Assuming that each component has two states of failure and success, and components failures are independent of each other. If the random value is less than the failure rate, the component will be unavailable. This simple way of modeling considers multiple-outage modes and is almost sufficient in a non-sequential simulation. When a non-sever failure occurs, some components can still be operated in different derated states. These states are not modeled in this work.

The probability of each sample can be calculated using Eq. (17). Finally the output of this step is a matrix, each column of which represents a sample of uncertain variable states.

(17) \mathrm{Prob}\left(sm\right)\mathrm{=Rand}\left(sm,r\right).

Scenario reduction

The computational requirements for solving scenario-based optimization models are directly affected by the number of scenarios. Therefore, an effective scenario reduction technique could be very essential and useful in solving large-scale problems. Reference [ 18] defines the reduction technique as a scenario-based approximation with a smaller number of scenarios and a reasonably good approximation of the original system. The scenario reduction technique that is applied in this study controls the goodness-of-fit of approximation by measuring a distance of probability distributions as a probability metric. After performing scenario reduction, a subset of scenarios with regarding probabilities is selected that models the initial probability distribution in terms of probability metrics. Efficient algorithms based on backward and fast forward methods are derived that determine optimal reduced measures. An overview of the simultaneous backward reduction method, based on Ref. [ 19], is given as follows.

Let ξ_s (s = 1,2,⋯,Sc) denote Sc different scenarios, each with a probability of Prob_s, and Dist_s,s′ be the distance of scenario pair (s, s′€). The simultaneous backward and fast forward is given in the following steps:

Step 1 Set S as the initial set of scenarios; DS is the set of scenarios to be omitted. The initial DS is null. Compute the distances of all scenario pairs:

Dist_s,s′ = Dist(ξ_s, ξ_s′), s, s′=1, 2, …, Sc;

Step 2 For each scenario k, Dist_k(r) = min Dist_k,s′, while s′, k\in S and s′ ≠ k, r is the index of scenario that has the minimum distance with scenario k;

Step 3 Compute DD_k(r) = Prob_k × Dist_k(r), k\in S . Choose d so that DD_d = min (DD_k), k\in S ; S = S-{d}, DS= DS+ {d}; Prob_r= Prob_r+ Prob_d;

Step 4 Repeat Steps 2–3 until the number to be deleted meets the predefined number of scenarios.

Proposed method

The aim of operation and planning in deregulated power systems is to maximize the social welfare through minimization of costs of the network, while the electric power is delivered to the customers with sufficient quality and reliability. Because of the high investment in DERs, there is considerable risk in their application. Therefore, the optimal placement and sizing of DERs are the most important steps to be performed considering various aspects of distribution networks. The objectives of this study are loss minimization, reduction of power which should be purchased from electricity market, loss reduction at the peak load level and improvement of voltage profile of the power system through proper application of DERs.

Objective fuzzification

Each objective in the fuzzy domain is associated with a membership function. The membership function specifies the degree of satisfaction of the objective. In the crisp domain, the objective is either satisfied or violated, indicating membership values of unity and zero, respectively. On the contrary, fuzzy sets consider varying degrees of membership function values from zero to unity [ 16]. The present work considers the following objectives for the DER placement problem.

1) Maximization of the saving by minimization of the energy loss, power purchased and loss at the peak load level due to the application of DERs.

2) Minimization of the voltage deviation at network buses.

Before continuing further in this section, the stochastic long-term load model should be reviewed to find out how to use it for DER placement.

In each scenario, there are 12 load levels representing four different load blocks of the three years of study horizon, each with a chance of occurrence. Each scenario itself has a probability. Combining the load levels of the scenarios, the total load probability density function (PDF) is obtained. The load PDF is divided into equal sections each with an occurrence probability. The centers of these sections are introduced into the optimization algorithm as the load levels.

The membership function consists of a lower and upper bound value along with a strictly monotonically decreasing and continuous function as described in Section 4.2.

Membership function for net saving

The net saving at kth load level due to application of DER in a distribution system is given in Eq. (18) (it should be noted that load levels are, in fact, the stochastic load levels that, along with line outages, reflect the system states in each scenario).

(18) N_S={\displaystyle \sum _{yr=\mathrm{1}}^{N_{yr}}{K}_{\mathrm{P}}{T}^{\mathrm{peak}}{\mathrm{LR}}_{yr}^{\mathrm{peak}}}-{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{DER}}}{K}_{\mathrm{inv}}^{\mathrm{DER}}{P}_i^{\mathrm{DER,Max}}}+{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{k}}}\left[K_{\mathrm{E}}{\rho }_k{\mathrm{LR}}_k+K_{\mathrm{E}}{\rho }_k{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{DER}}}{P}_{i,k}^{\mathrm{DER}}}-{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{DER}}}{K}_k^{\mathrm{DER}}{P}_{i,k}^{\mathrm{DER}}}\right]},

where K_P is a factor to convert peak power loss reduction to dollar ($/kW); T^peak is the duration of peak load (hours); {\mathrm{LR}}_{yr}^{\mathrm{peak}} is the power loss reduction at peak load level at year yr due to application of DERs (kW); K_E is a factor to convert energy losses to dollar ($/kWh); {\rho }_k is the probability of kth load level; LR_k is the reduction in power loss at kth load level due to application of DERs (kWh); N_DER is the number of DERs; P_{i,k}^{\mathrm{DER}} is the power output of ith DER at kth load level (kWh); K_k^{\mathrm{DER}} is the cost of operation and maintenance of DER at kth load level ($/kWh); K_{\mathrm{inv}}^{\mathrm{DER}} is the cost of investment of DER ($/kW); and P_i^{\mathrm{DER,Max}} is the maximum capacity of ith DER (kW).

Considering a positive profit for application of DERs, for net saving in (13), there is N_S>0 which means

(19) N_{\mathrm{S}}={\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{k}}}\left[K_{\mathrm{E}}{\rho }_k{\mathrm{LR}}_k+K_{\mathrm{E}}{\rho }_k{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{DER}}}{P}_{i,k}^{\mathrm{DER}}}\right]}+{\displaystyle \sum _{yr=\mathrm{1}}^{N_{yr}}{K}_{\mathrm{P}}{T}^{\mathrm{peak}}{\mathrm{LR}}_{yr}^{\mathrm{peak}}-{\displaystyle \sum _{k=\mathrm{1}}^{N_k}{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{DER}}}{K}_k^{\mathrm{DER}}{P}_{i,k}^{\mathrm{DER}}}}}\ge \mathrm{0},

(20) {\displaystyle \frac{{\displaystyle \sum _{k=\mathrm{1}}^{N_k}{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{DER}}}{K}_k^{\mathrm{DER}}{P}_{i,k}^{\mathrm{DER}}}}}{{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{k}}}\left(K_{\mathrm{E}}{\rho }_k{\mathrm{LR}}_k+K_{\mathrm{E}}{\rho }_k{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{DER}}}{P}_{i,k}^{\mathrm{DER}}}\right)}+{\displaystyle \sum _{yr=\mathrm{1}}^{N_{yr}}{K}_{\mathrm{P}}{T}^{\mathrm{peak}}{\mathrm{LR}}_{yr}^{\mathrm{peak}}}}}\le \mathrm{1.}

Define

(21) x_k={\displaystyle \frac{{\displaystyle \sum _{k=\mathrm{1}}^{N_k}{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{DER}}}{K}_k^{\mathrm{DER}}{P}_{i,k}^{\mathrm{DER}}}}}{{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{k}}}\left(K_{\mathrm{E}}{\rho }_k{\mathrm{LR}}_k+K_{\mathrm{E}}{\rho }_k{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{DER}}}{P}_{i,k}^{\mathrm{DER}}}\right)}+{\displaystyle \sum _{yr=\mathrm{1}}^{N_{yr}}{K}_{\mathrm{P}}{T}^{\mathrm{peak}}{\mathrm{LR}}_{yr}^{\mathrm{peak}}}}}.

Equation (15) indicates that if x_k is high, the net saving (profit) is low and vice versa. The membership function for net saving (profit) is given in Fig. 3. Based on Fig. 3, Eq. (22) can be obtained.

(22) \mu s_k={\displaystyle \frac{x_{\max }-x_k}{x_{\max }-x_{\min }}}\mathrm{for} x_{\min }\le x_k\le x_{\max },

\mu s_k=\mathrm{1}\mathrm{for} x_k\le x_{\min },

(23) \mu s_k=\mathrm{0}\mathrm{for}{x}_k\ge x_{\max }.

In this study, x_max is assumed to be 1.0; in order to achieve x_min the proposed method is once performed without consideration of voltage improvement as one of the objectives. The value of x_min is determined based on the maximum profit to cost ratio. This means that if the maximum profit to cost ratio achieved is 0.6, x_min will be 0.375. The x_min of 0.375 means that the unity membership value is assigned if the savings is 37.5% or more, and the x_max of 1.0 means zero membership value is assigned if the profit is zero percent of the cost or has a negative value.

Membership function for node voltage deviation

The basic purpose of this membership function is that the deviation of node voltage should be minimized. At kth load level of the load duration curve, define

(24) y_k=\max \left[\mathrm{abs}\left(V_{\mathrm{n}}-V_{i,k}\right)\right]\mathrm{for}i=2,3,\cdots ,N_{\mathrm{B}},

where V_{i,k} is the voltage magnitude of node i at kth load level per unit and V_{\mathrm{n}} is the nominal voltage magnitude that is equal to 1 per unit. It should be noted that it is for the case that the substation is located at bus 1. As the maximum value of nodes voltage deviation decreases, the assigned membership value increases and vice versa. Figure 4 illustrates the membership function for maximum node voltage deviation defined in Ref. [ 16]. Based on Fig. 3 and taking into account y_{\max }^{-}\le y_k\le y_{\max }^{+} , there are

(25) \mu V_k=\{\begin{array}{l}{\displaystyle \frac{y_{\max }^{+}-y_k}{y_{\max }^{+}-y_{\min }^{+}}}{y}_k\ge \mathrm{0}\hfill \\ {\displaystyle \frac{y_{\max }^{-}-y_k}{y_{\max }^{-}-y_{\min }^{-}}}{y}_k\le \mathrm{0}\hfill \end{array},

(26) \mu V_k=\mathrm{1}\mathrm{for}{y}_{\min }^{-}\le y_k\le y_{\min }^{+},

(27) \mu V_k=\mathrm{0}\mathrm{for}{y}_k\ge y_{\max }^{+}\mathrm{or}{y}_k\le y_{\max }^{-}.

In this study, y_{\min }^{+} and y_{\max }^{+} are considered to be 0.05 and 0.10, respectively. y_{\min }^{-} and y_{\max }^{-} are assumed to be –0.05 and –0.10, respectively. Considering V_n = 1, y_{\min }^{+} = 0.05 means the minimum system voltage will be 0.95 pu, which means that if the minimum system voltage is greater than or equal to 0.95 pu, the membership value is 1. Similarly, y_{\max }^{+} = 0.10 means the minimum system voltage allowed will be 0.90 pu and if the minimum system voltage is less than or equal to 0.90 pu, the assigned membership value will be 0. For more information about practical issues concerning the choice of the membership functions and weighting factors, please refer to Refs. [ 12– 14].

Fuzzy formulation for several objectives

The two fuzzified objectives described in the previous section are dealt with by integrating them into a fuzzy satisfaction objective function F, through appropriate weighting factors (K) as Eq. (23).

(28) \mathrm{Max}F=K_{\mathrm{1}}\mu s_k+K_{\mathrm{2}}\mu V_k,

where K₁ and K₂ are weighting factors considered in this study for investigation of the impact of each objectives in planning of the DERs. K_{\mathrm{1}} = 0 and K_{\mathrm{2}} = 1 means that only the improvement of voltage profile is considered as the objective of optimization and vice versa; while K_{\mathrm{1}} = 0.5 and K_{\mathrm{2}} = 0.5 means that these two objectives are assumed to be equally important. The weighting factors can be determined according to the preferences of the operators [ 16].

Simulation results

To test the effectiveness of the proposed method, the IEEE 30-bus distribution system including 22 load points, 6 auxiliary substations and a main feeder [ 20] is chosen. This system is shown in Fig. 5.

The load pattern in the peak load level of the present year is listed in Table 1. Both real (kW) and reactive (kvar) loads are specified. The base value of voltage and power are 23 kV and 100 MVA, respectively. Table 1 also lists the resistance and reactance of the lines. After scenario aggregation, ten load levels are considered, each with an occurrence probability. These load levels are presented in Table 2 as the percentage of the peak load level of the present year. More load levels can be considered for the sake of accuracy. The IHBMO parameters are tabulated in Table 3.

It is assumed that the size of DERs varies in 100 kW steps. The investment and operation costs of DERs are borrowed from Ref. [ 21].

Fuzzy optimization problem considering several objectives, deterministic case

Before testing the proposed stochastic method, a deterministic version of the proposed method is tested on the IEEE 30-bus distribution system in this sub-section. The solution to this deterministic problem can be compared with the solution to the stochastic problem to show the necessity of the stochastic modeling of the problem. The values of the load and energy growth rates are considered to be 0.08. The load model with nine levels is used for the three-year time horizon. Table 4 shows these load levels and the regarding time durations. The proposed stochastic approach can be simply modified for the deterministic problem by substituting the probability of load levels ( {\rho }_k ) with load level durations in Table 4.

To find the maximum attainable profit which, as discussed earlier, is an important factor in construction of membership functions, first, the voltage deviation is omitted from the objective function and the maximum profit found is 1389923 $. Now, the suitable membership function regarding to profits can be found to solve the problem. Table 5 shows the optimal solution of the deterministic problem.

Figure 6 displays the voltage profile for compensated and uncompensated systems for the peak load level of the present year. As can be seen from Fig. 6 that the voltage deviation is lower for the compensated system, which demonstrates that the algorithm can effectively mitigate the voltage deviation while the value of profit is still acceptable compared with the maximum attainable profit gained in previous case study.

Fuzzy optimization problem considering several objectives, stochastic case

In this case study, the proposed stochastic approach is used to find the best solution to the optimization problem considering several objectives. To find the shape of membership function of the first part of the objective function, initially the stochastic problem is solved considering the profit as the objective function to find the maximum attainable value of profit. Table 6 shows the results of the placement problem for the single objective problem. As can be seen from Table 6, the maximum attainable profit with the minimum acceptable BCR is 6658313.1 $. At the next stage the stochastic problem is solved considering all the objectives and the results, as given in Table 7.

Figure 7 demonstrates the convergence characteristic for the stochastic problem by IHBMO. It depicts the change in the BCR of the best solution (Queen) versus iterations of the algorithm. As it can be seen, the IHBMO has rapid convergence characteristic.

Statistical analysis of the results

Since the algorithm used in this paper is a heuristic optimization algorithm and the stochastic nature of the system has been taken into account, the results derived from the proposed method might vary in each run. To investigate the effect of these factors, statistical analysis of the results is discussed in this subsection. The proposed algorithm is run 50 times to find the best solution of the problem discussed in case study B (stochastic case with both objectives included).

Table 8 shows the statistical analysis of the results. The IHBMO parameters are considered to be fixed in all runs. As can be seen, the standard deviation of the solutions is very low. This shows the robustness of the proposed algorithm against some factors such as initial population of IHBMO algorithm. Another study is also conducted to analyze the effects of the increase in degree of uncertainty associated with system loads. The standard deviation of peak load and energy growth is changed from 2% to 5% and the results are presented in Table 9. As can be seen from Table 9, the profit decreases as the degree of uncertainty increases. It is also interesting to observe that as the degree of uncertainty increases, the maximum value of voltage deviation increases with one exception. The reason lay under this fact may be that, with the increase in these standard deviations, the total objective function definitely decreases, but each objective may show unexpected trend.

Discussion

As said earlier, the proposed algorithm schedules the output of DERs in each load levels individually to avoid the unwanted rejection of optimal solutions. To illustrate this, consider the result of Table 7. The reduction in the power purchased (PPR) from the electricity market simply shows that in all load levels the maximum output of DERs is not scheduled. The maximum power which can be supplied by DERs is 7900 kW, so the maximum energy which can be supplied is 3 × 8760 × 7.9= 207612 MWh while the energy served by DERs is 197015 MWh. Again the minimum BCR of 1.3 is considered as a constraint.

Comparing the results of Table 7 with those presented in Table 6, the profit and the maximum voltage deviation are less for the problem in which several objectives are considered. This shows that the improvement of voltage profile causes a small reduction in profit.

Comparing the results of stochastic algorithm (Table 7) with those obtained in the case study (Sub section 5.1) for the deterministic problem, the necessity of the stochastic modeling of the problem can be understood. To clarify this point, the solution of deterministic problem is fed into the stochastic model to calculate the stochastic profit gained by applying the best solution to the deterministic problem. The profit value is 4412295.2 $ which is so much lower than the stochastic profit gained by best solution of stochastic problem reported in Table 7 (6255752.6 $).

The voltage profile in the peak load level of each year of the three-year time horizon of the study for compensated and uncompensated cases are presented in Fig. 8. As can be seen from Fig. 8, the proposed method can effectively mitigate the voltage deviation. It should be noted that the optimal power flow algorithm for the peak load level of years 2 and 3 in uncompensated case does not converged with the main feeder voltage of 1.05 pu; so the voltage of main feeder for these states is considered to be 1.1 pu. Table 10 shows the mean deviation of voltage at different load points and in fact summarized the results presented in Fig. 8.

Conclusions

An IHBMO based fuzzy stochastic long-term optimization methodology considering several objectives has been proposed in this paper for optimal placement and sizing of DERs. As the results of case studies shows, ignoring the uncertainties in DER placement problem renders the stochastic saving gained non-optimal. The optimization algorithm simultaneously seeks the reduction in power loss; power purchased from the electricity market, power loss at the peak load level and deviation of the voltage magnitude at the load points. A proper modeling of economic aspects of the problem has also been presented in this paper. The results of case studies show that the proposed algorithm can effectively guarantee the justification of investment on DERs. The results also indicate that the method schedules the output power of DERs in each load level and avoids the inconvenience rejection of more optimal solutions. In the fuzzifying process, the membership functions have been obtained with an effective method, instead of the predefined membership functions. Some other aspects of the proposed method have also been discussed in case studies.

References

[1] Pai P F, Yang S L, Chang P T. Forecasting output of integrated circuit industry by support vector regression models with marriage honey-bees optimization algorithms. Expert Systems with Applications, 2009, 36(7): 10746–10751

[2] Ghasemi A. A fuzzified multi objective interactive honey bee mating optimization for environmental/economic power dispatch with valve point effect. International Journal of Electrical Power and Energy Systems, 2013, 49: 308–321

[3] Shayeghi H, Ghasemi A.Multiple PSS design using an improved honey bee mating optimization algorithm to enhance low frequency oscillations. International Review of Electrical Engineering (I.R.E.E.), 2011, 6(7): 3122–3133

[4] Jain N, Singh S N, Srivastava S C. A generalized approach for DG planning and viability analysis under market scenario. IEEE Transactions on Industrial Electronics, 2013, 60(11): 5075–5085

[5] Gil H A, Joos G. Models for quantifying the economic benefits of distributed generation. IEEE Transactions on Power Systems, 2008, 23(2): 327–335

[6] Chiradeja P, Ramakumar R. An approach to quantify the technical benefits of distributed generation. IEEE Transactions on Power Systems, 2004, 19(4): 764–773

[7] Prenc R, Škrlec D, Komen V. Distributed generation allocation based on average daily load and power production curves. International Journal of Electrical Power & Energy Systems, 2013, 53, 612–622

[8] Wang C, Nehrir M H. Analytical approaches for optimal placement of distributed generation sources in power systems. IEEE Transactions on Power Systems, 2004, 19(4): 2068–2076

[9] Georgilakis P S, Hatziargyriou N D. Optimal distributed generation placement in power distribution networks: models, methods, and future research. IEEE Transactions on Power Systems, 2013, 28(3): 3420–3428

[10] El-Khattan W, Bhattacharya K, Hegazy Y, Salama M M A. Optimal investment planning for distributed generation in a competitive electricity market. IEEE Transactions on Power Systems, 2005, 20(4): 1718–1727

[11] Ahmadigorji M, Abbaspour A, Rajabi-Ghahnavieh A, Fotuhi-Firuzabad M. Optimal DG placement in distribution systems using cost/worth analysis. World Academy of Science, Engineering and Technology, 2009, 25: 746–753

[12] Wang S, Watada J, Pedrycz W. Recourse-based facility location problems in hybrid uncertain environment. IEEE Transactions on Systems, Man, and Cybernetics. Part B, Cybernetics, 2010, 40(4): 1176–1187

[13] Wang S, Watada J, Pedrycz W. Value-at-Risk-based two-stage fuzzy facility location problems. IEEE Transactions on Industrial Informatics, 2009, 5(4): 465–482

[14] Wang S, Watada J. A hybrid modified PSO approach to VaR-based facility location problems with variable capacity in fuzzy random uncertainty. Information Sciences, 2012, 192(1): 3–18

[15] Celli G, Ghiani E, Mocci S, Pilo F. A multiobjective evolutionary algorithm for the sizing and sitting of distributed generation. IEEE Transactions on Power Systems, 2005, 20(2): 750–757

[16] Das D. Optimal placement of capacitors in radial distribution system using a Fuzzy-GA method. International Journal of Electrical Power and Energy Systems, 2008, 30(6–7): 361–367

[17] Neimane V. Distribution network planning based on statistical load modeling applying genetic algorithms and Monte-Carlo simulations. In: IEEE Power Tech Proceedings. Porto, Portugal, 2001

[18] Roh J H, Shahidehpour M, Fu Y. Market-based coordination of transmission and generation capacity planning. IEEE Transactions on Power Systems, 2007, 22(4): 1406–1419

[19] Hashemzadeh H, Ehsan M. Locating and parameters setting of unified power flow controller for congestion management and improving the voltage profile. Asia-Pacific Power and Energy Engineering Conference. Chengdu, China, 2010

[20] Santoso N I, Tan O T. Neural-net based real-time control of capacitors installed on distribution system. IEEE Transactions on Power Delivery, 1990, 5(1): 266–272

[21] Basu A K, Chowdhury S, Chowdhury S P. Impact of strategic deployment of CHP-based DERs on microgrid reliability. IEEE Transactions on Power Delivery, 2010, 25(3): 1697–1705

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg