# A Sustainability-Oriented Multiobjective Optimization Model for Siting and Sizing Distributed Generation Plants in Distribution Systems

Authors: Guang Chen; Bin Chen; Pan Dai; Hao Zhou
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper proposes a sustainability-oriented multiobjective optimization model for siting and sizing DG plants in distribution systems. Life cycle exergy (LCE) is used as a unified indicator of the entire system’s environmental sustainability, and it is optimized as an objective function in the model. Other two objective functions include economic cost and expected power loss. Chance constraints are used to control the operation risks caused by the uncertain power loads and renewable energies. A semilinearized simulation method is proposed and combined with the Latin hypercube sampling (LHS) method to improve the efficiency of probabilistic load flow (PLF) analysis which is repeatedly performed to verify the chance constraints. A numerical study based on the modified IEEE 33-node system is performed to verify the proposed method. Numerical results show that the proposed semilinearized simulation method reduces about 93.3% of the calculation time of PLF analysis and guarantees satisfying accuracy. The results also indicate that benefits for environmental sustainability of using DG plants can be effectively reflected by the proposed model which helps the planner to make rational decision towards sustainable development of the distribution system.


---

## Body

## 1. Introduction

Distributed generation (DG) is developing fast all over the world in recent years due to its promising potential to reduce the portion of fossil energy consumption in electric power generation and mitigate power losses and harmful emissions [1, 2]. For the distribution system, DG plants can provide valuable reserve capacities and improve the flexibility of system operation [3–5]. However, if not properly planned, integrating DG plants may also cause some negative impacts, especially when they are driven by renewable energies such as wind and solar power, since these renewable energies are generally unstable and intermittent [6]. Additional difficulties in reactive power compensation, load flow control, fault protection, and power quality regulation are other possible negative impacts [7].The performance of the DG plants in a distribution system is greatly impacted by the arrangement, that is, the sites and sizes, of the plants in the system [1, 8]. Thus, it is necessary to optimize the arrangement in the planning stage since it is much more difficult to modify the arrangement thereafter. Such optimization process is also called siting and sizing of DG plants [9]. Presently, several approaches have been proposed for siting and sizing DG plants in distribution systems, and brief reviews can be found in [1, 10]. The present optimization approaches can be categorized into analytical methods [11–13], metaheuristic methods [14, 15], artificial intelligent methods [16–19], and others [20, 21]. The technical and economic issues, for example, power losses, node voltages, line capacities, and economic costs, are mainly concerned in previous studies. Comparatively, the environmental sustainability, including various environmental effects such as harmful emissions and natural resource depletion of using DG plants in a distribution system is not fully evaluated. Most studies tend to presume the effectiveness and equivalence of different DG plants in environment protection and to maximize the installed capacities considering some technical and economic constraints. However, the environmental impacts of different DG plants may vary considerably due to the differences in primary energy sources, production techniques, and utilization efficiencies of the generators. Thus, it is necessary and urgent to properly evaluate the environmental effects in the siting and sizing process to obtain better sustainability. Life cycle assessment (LCA) is generally regarded as an excellent measure to evaluate environmental effects of an energy system during its life cycle [22], but it lacks a mechanism to comprehensively measure various types of environmental factors. The weighting aggregation and the so-called neoclassical economic theory are two typical options for compensation, but there are some drawbacks of these methods, such as subjectivity, incompleteness, and overoptimism [23]. Exergy analysis is another option that takes exergy as a general indicator to evaluate the various factors and has the potential to overcome the above drawbacks [22]. Moreover, its unique capability of identifying the irreversibilities of energy systems by the second law of thermodynamics makes exergy an excellent indicator to measure the system’s long-term sustainability [24, 25]. Several studies have combined exergy analysis with LCA to analyze or optimize the designs of energy systems [22, 26–28], forming the so-called exergetic LCA (ELCA). In this paper, the ELCA method is used to build a unified quantitative indicator of environmental sustainability for installing and utilizing DG plants in a distribution system, and then the siting and sizing problem is modeled as a multiobjective chance-constrained integer optimization problem. Thus, the proposed optimization model not only minimizes the economic cost and power loss but also optimizes environmental sustainability.Another key problem in the optimization process is how to efficiently reflect the performance of the distribution system considering the unstable outputs of renewable energies, load levels, and other uncertain factors [7]. Presently, probabilistic load flow (PLF) analysis is widely used due to its capability of reflecting static system security [1, 7, 10, 18, 19]. Improving the efficiency of PLF analysis is of great significance since it is performed repeatedly for every candidate solution to determine its feasibility and fitness in the optimization process. Studies have proposed several methods for PLF analysis, such as cumulant-based series expansion methods (including Gram-Charlier expansion, Edgeworth expansion, and Cornish-Fisher expansion) [29], point estimate method [30], and sampling-based simulation methods (e.g., Monte Carlo sampling (MCS) [29, 30], Latin hypercube sampling (LHS) [31], etc.). In general, the sampling-based simulation methods are more accurate and reliable than the cumulant-based series expansion methods and the point estimate method but need much more calculation time [29–31]. Compared with MCS which samples randomly, LHS uses a more complicated sampling mechanism to obtain semirandom or deterministic samples and thus improves efficiency while guaranteeing accuracy [31–33]. To further improve the efficiency of LHS-based PLF analysis, a semilinearized simulation method is proposed to combine with it in this paper. In the proposed method, the time-consuming iteration process for nodal voltage calculation is avoided by using  linearized  load flow equation. Satisfying accuracy and reliability is guaranteed by retaining the nonlinearity of power flow calculation and the framework of LHS-based PLF analysis.The rest of this paper is organized as follows: the models and methods are described in Section2; Section 3 gives a numerical case study based on the modified IEEE (Institute of Electrical and Electronics Engineers) 33-node system; Section 4 concludes.


## 2. Methodology

### 2.1. Exergetic Life Cycle Assessment

ELCA evaluates all the exergies consumed within the life cycle of a system. As shown in (1), the life-cycle exergy (LCE) of a DG plant can be calculated as the sum of the cumulative exergy consumption (CEC) and the emission abatement exergy (AbatE). Here, CEC is defined as the sum of the exergies of all the consumed resources and reflects resource depletion [22]. AbatE is defined as the cumulative equivalent exergy caused by primary resource consumption to remove or separate the wasted emissions by an implemented or proposed treatment process and reflects the environmental effect of emissions [23, 34, 35]. Thus, the environmental sustainability of using a DG plant can be reflected by LCE [22–24, 27]:
(1)LCE=CEC+AbatE.For convenience, the life cycle of a DG plant can be divided into 3 phases, that is, construction phase, operation phase, and disposal/recycle phase. Then, CEC and AbatE can be calculated as(2)CEC=CECCon+CECOpr+CECDis,AbatE=AbatECon+AbatEOpr+AbatEDis,
where Con, Opr, and Dis represent the construction, operation, and disposal/recycle phases, respectively.In the construction phase, the cumulative exergy consumption (CECCon) is mainly caused by the consumption of fuels and the inputs of nonfuel materials (e.g. iron ore, limestone, sand, etc.) needed for constructing the DG plant. The exergies of fuels and non-fuel materials can be calculated by multiplying the consumed/input amounts with corresponding exergy content factors, that is,
(3)CECCon=CDG·(∑iαiqiaCon,i+∑jxjbCon,j),(4)AbatECon=CDG∑kβkeCon,k,
where CDG is the installed capacity of the DG plant; αi is the ratio between exergy and energy content for fuel i which can be collected from relevant studies; qi is the energy content of per unit of fuel i; aCon,i is the amount of fuel j consumed to construct per capacity of the DG plant; xj is the chemical exergy content of the non-fuel material j; bCon,j is the mass of the non-fuel material j needed for per capacity of the DG plant; eCon,k is the emission amounts of pollutant k in the construction stage and the operation stage, respectively; βk is the abatement exergy factor of pollutant k.For the operation phase, the input of non-fuel materials is negligible, and thus the cumulative exergy consumption is mainly caused by fuel consumption which is proportional to power generation, shown as (5). For the DG plants that solely use renewable energies, the CECOpr can be neglected. Similarly, the emission abatement exergy of fuel-driven DG can be calculated by multiplying the amounts of different emissions and abatement exergy factors as shown in (6) and that of the renewable energy-driven DG is neglected,
(5)CECOpr=CDGHDGTDG∑iαiqiaOpr,iηOpr,i,(6)AbatEOpr=CDGHDGTDG∑kβkeOpr,k,
where HDG and TDG are the annual utilization hours (AUH) and the operation expectancy of the DG plant, respectively; aOpr,i is the cumulative mass of fuel j consumed to generate per unit of electric power; ηOpr,i is the overall efficiency of production and delivery for fuel i. eOpr,k is the emission amount of pollutant k in the construction stage and the operation stage, respectively.The cumulative exergy consumption and emission abatement exergy in the disposal/recycle phase consist of exergy consumption for unit disposal and exergy saving for unit recycle. In general, these exergies can be calculated proportionally to the exergies in the construction phase, shown as follows:(7)CECDis=-λCECCECCon,AbatEDis=-λAbatAbatECon,
where λCEC and λAbat represent the exergy saving ratios with respect to CECCon and AbatECon, respectively.


### 2.2. Optimization Model

The optimization variables of the siting and sizing problem are the locations of the DG plants in the studied distribution network and their installed capacities. The location variables can be expressed with integer numbers. Usually a DG unit is produced with a standard capacity, and thus the capacity variables of a DG plant can be expressed as integers that represent the number of installed DG units. Based on such arguments, the siting and sizing problem is modeled as a multi-objective integer programming problem in this paper. The objectives include minimizing the economic cost, life-cycle exergy consumption, and power losses.


#### 2.2.1. Objective Functions

The multi-objective function can be written as(8)MinF={F1,F2,F3},
where F1 represents the life-cycle economic cost of utilizing DG plants; F2 represents the life-cycle exergy of using DG plants; F3 represents the expected power loss considering the uncertainties of power loads and DG outputs.The life-cycle economic cost of utilizing DG plants can be calculated as(9)F1=∑i∑j(CCon,ij+COpr,ij+CDis,ij),
where CCon,ij, COpr,ij and CDis,ij, are economic costs spent on construction, operation, and disposal/recycle to use jth type of DG at node i which can be calculated as
(10)CCon,ij=cCon,jCDG,ij,COpr,ij=ξCPV,jcOpr,jCDG,ijHDG,j,CDis,i=-ζTDG,j+1dCost,jCCon,ij,
where cCon,j is the economic cost for constructing per capacity of jth-type DG; cOpr,j is the economic cost for generating per unit of electricity by jth-type DG; dCost,j is the recycle coefficient for jth-type DG; coefficients ξCPV,j and ζ can be calculated as (11), in which f and δ stand for the inflation rate and discount rate, respectively. Generally, due to the intermittency of renewable energies, both the operation costs and power generation of DG plants that solely use renewable energies are much less than those of the DG plants which have the same capacities but use fuels:
(11)ξCPV,j=ζ-ζTDG,j+11-ζ,ζ=1+f1+δ.The objective function of minimizing life-cycle exergy can be written as(12)F2=∑i∑j(LCEDG,ij-ENoDG,ij),
where LCEDG,ij is the life-cycle exergy used by the DG plant of jth type at node i; ENoDG,ij stands for the expected exergy consumption to provide the same amount of DG-generated energies (electricity, heat, refrigeration, etc.) by alternative sources if DG plant ij is not installed.E
N
o
DG
,
i
j can be calculated as
(13)ENoDG,ij=CDG,ijHDG,j∑t=1TDG,j(γE,t+∑kυjkγk,t),
where γE,t represents the exergy needed for providing per unit of electric power in tth year by alternative sources; υjk is the amount of by-product k when generating per unit of electric power; γk,t stands for the exergy needed for providing per unit of by-product k in tth year by alternative sources.The objective function of minimizing expected power loss as(14)F3=E{ELoss},
where ELoss stands for the power loss and E{ELoss} represents the mean value of ELoss.


#### 2.2.2. Constraints

The total capacity of installed DG plants should be no higher than a given cap for the sake of safety and reliability of power supply, that is,(15)∑i∑jCDG,ij≤CDGmax.Other main constraints include the nodal voltage magnitude limitations and the branch capacity limitations. In real operation, these limitations can be violated for short terms under some operating conditions. Considering the uncertainties of renewable energy sources and power loads, chance constraints are generally used to guarantee the practicability of the solution shown as follows:(16)Pr{Vimin<Vi<Vimax}≥ηV,Pr{|Si,j|<Si,jmax}≥ηS,
where Si,j stands for the apparent power flow in the branch from node i to node k; Vimin, Vimax, and Si,jmax are limitation values; ηV and ηS are required confidence parameters.For each operation condition, the load flow equations need to be obeyed, that is,(17)Pi=Vi∑j∈iVj[Gijcos(θi-θj)+Bijsin(θi-θj)],Qi=Vi∑j∈iVj[Gijsin(θi-θj)-Bijcos(θi-θj)],
where j∈i represents the set of every node j which is directly connected with node i by a branch in the distribution network, including node i itself; Pi and Qi stand for the active and reactive powers injected into node i; Vi and θi stand for the voltage magnitude and angle of node i; Gij and Bij represent the real and imaginary parts of admittance between node i and j.


### 2.3. Treatment of Uncertainties

#### 2.3.1. Models of Uncertain Variables

Wind Turbine Output. Generally, the stochastic behaviors of the wind speed v in a wind farm can be described by the Weibull distribution whose probability density function (PDF) can be written as [37]
(18)fv(v)=kckvk-1exp[-(vc)k],
where k and c are shape factor and scale factor of the wind speed PDF at the wind farm.Ignoring minor nonlinearities, a turbine’s power output can be calculated with a given wind speedv:
(19)Pw={0v<vIorv>vOPwRv-vIvR-vIvI≤v≤vRPwRvR<v≤vO,
where vI, vO and vR are cut-in, cut-out, and rated wind speeds for the wind turbine; PwR is the turbine’s rated output.PV Output. Ignoring the minor factors, the output of a PV unit can be characterized as a function of solar irradiance and PV cell temperature [29, 36]:
(20)Ps=rrmax(1-κ·ΔT)Ps,max,
where r and rmax stand for the solar irradiance and its maximum value; Ps,max represents the maximum output at the rated cell temperature; ΔT is the deviation of the PV cell’s temperature from the rated one; κ is a coefficient.In the daytime, the PDF of solar irradiance can be characterized by Beta distribution [18, 29] which can be written as
(21)fr(r)=Γ(αr+βr)Γ(αr)Γ(βr)(rrmax)αr-1(1-rrmax)βr-1,
where αr and βr are distribution coefficients.The PV cell temperature is determined by the ambient temperature and prevailing wind speed. The deviation of the PV cell’s temperature can be described by normal distribution [29], that is,
(22)ΔT~N(0,σΔT2).Power Load. Generally, the uncertainty of the power load at a node can be described with normal distribution while the power factor can be regarded as a constant for simplicity [7, 38]. Therefore,
(23)PL,i~N(μL,i,σL,i2),
where PL,i represents the active power load at node i; μL,i and σL,i2 are the mean value and variance of PL,i.LHS-Based Semilinearized Probabilistic Load Flow Analysis. The typical procedure of LHS-based probabilistic load flow analysis can be found in [31], which consists of 3 steps: sampling, simulating, and accounting. In the sampling step, a required number of sample arrays of the uncertain variables are generated according to the variables’ probability density functions (see Section 2.3.1). Correlations of variables can be taken into account by Nataf transformation [31]. Each sample array represents a possible operating point and is used to evaluate the corresponding node voltages and power flows by solving (17) in the simulating step. Accounting on simulation results is then performed to obtain the statistical properties of certain variables, such as the probabilities in (16) and the expected power loss in (14).For simplicity, the load flow equation, that is, (17), can be expressed by (24) in which vector Y expresses the injected powers and vector X stands for the nodal voltages. Then, the power flows in branches can be expressed as a function of X.  Given an initial operating point Y0, the corresponding X0 and S0 can be calculated. However, since time-consuming iterations are needed due to the complexity and nonlinearity of (24), that is, (17), calculating X0 needs much more time than calculating S0 by (25), and this makes the simulating step requires much more time than the other 2 steps in PLF analysis:
(24)Y=g(X),(25)S=h(X).Expanding (24) at Y0 according to Taylor series and omitting the terms higher than the second order will yield
(26)X=X0+ΔX=X0+M0ΔY,
where vector ΔY expresses the changes of Y from Y0; vector ΔX stands for the changes of X; M0=J0-1, and J0 is the final Jacobian matrix at the initial operating point Y0. Equation (26) is also called “linearized load flow equation.”Similarly, (25) can be linearized as
(27)S=S0+ΔS=S0+T0ΔY,
where T0=L0M0 and L0=(∂h(X))/∂X|X=X0.Apparently, load flow simulation based on (26) and (27) will need much less calculating time than based on (24) and (25), and this is the theoretical fundament of cumulant-based series expansion methods for PLF analysis [29]. However, simulations show that (26) generates quite accurate results even with large values of ΔY; while, on the contrary, (27) can not always produce accurate results since it needs two times linearization from ΔY to ΔX. Nevertheless, (26) can be utilized to improve the calculation efficiency of the simulating step in LHS-based PLF analysis following a semi-linearized simulation method as below.Step 1.
Select the operating condition when all the uncertain variables are fixed at their mean values as the initial operating pointY0. Substitute Y0 into (24) and solve the corresponding load flow equations to get the initial nodal voltage vector X0 and the Jacobian matrix J0.Step 2.
Seti=1, and calculate ΔYi=Yi-Y0, where Yi is the ith sample array generated by LHS method in the sampling step.Step 3.
SubstituteΔYi into (26) and calculate Xi.Step 4.
SubstituteXi into (25) to calculate Si.Step 5.
Seti=i+1, repeat Step 2 to Step 5 until all the sample arrays are evaluated.Thus, the semi-linearized simulation method avoids the time-consuming iterations in calculating the nodal voltages for different operating condition samples and guarantees the accuracy by keeping the nonlinearity of (25). Meanwhile, by using LHS-based PLF analysis framework, the inaccuracy and instability of cumulant-based series expansion are avoided in the PLF analysis.


### 2.4. Solving Approach and Methodology Framework

To solve the proposed multi-objective optimization model, the nondominated sorting genetic algorithm II (NSGA-II) [39] is used in this paper. NSGA-II greatly improves the performances of NSGA by introducing elitism and diversity-preserving mechanism algorithms and assures better convergence without losing solution diversity [40]. It should be noted that, instead of giving one optimal solution, NSGA and NSGA-II generate a set of nondominated solutions, that is, Pareto-optimal front, for multi-objective optimization problems according to the Pareto optimality principle [41]. Each solution actually represents the optimal solution under a different tradeoff relationship of the objectives. This enables the planner to see different possible solutions and helps the planner to choose a final solution based on certain tradeoff criteria according to a specific situation after comparing possible alternatives. The brief framework of the optimization algorithm used in this paper is shown in Figure 1, where the step of evaluating chromosomes consists of several substeps shown in Figure 2.Figure 1
Brief framework of the NSGA-II algorithm in this paper.Figure 2
Substeps in the step of evaluating chromosomes.Each chromosome in the algorithm represents a candidate solution for the optimization problem and is coded with integer genes as shown in Figure3, where each pane represents a gene in the chromosome. The position of the gene indicates the candidate site and type of a DG plant, and the value of the gene stands for the number of installed units.Figure 3
Coding scheme of the genes in a chromosome.


## 2.1. Exergetic Life Cycle Assessment

ELCA evaluates all the exergies consumed within the life cycle of a system. As shown in (1), the life-cycle exergy (LCE) of a DG plant can be calculated as the sum of the cumulative exergy consumption (CEC) and the emission abatement exergy (AbatE). Here, CEC is defined as the sum of the exergies of all the consumed resources and reflects resource depletion [22]. AbatE is defined as the cumulative equivalent exergy caused by primary resource consumption to remove or separate the wasted emissions by an implemented or proposed treatment process and reflects the environmental effect of emissions [23, 34, 35]. Thus, the environmental sustainability of using a DG plant can be reflected by LCE [22–24, 27]:
(1)LCE=CEC+AbatE.For convenience, the life cycle of a DG plant can be divided into 3 phases, that is, construction phase, operation phase, and disposal/recycle phase. Then, CEC and AbatE can be calculated as(2)CEC=CECCon+CECOpr+CECDis,AbatE=AbatECon+AbatEOpr+AbatEDis,
where Con, Opr, and Dis represent the construction, operation, and disposal/recycle phases, respectively.In the construction phase, the cumulative exergy consumption (CECCon) is mainly caused by the consumption of fuels and the inputs of nonfuel materials (e.g. iron ore, limestone, sand, etc.) needed for constructing the DG plant. The exergies of fuels and non-fuel materials can be calculated by multiplying the consumed/input amounts with corresponding exergy content factors, that is,
(3)CECCon=CDG·(∑iαiqiaCon,i+∑jxjbCon,j),(4)AbatECon=CDG∑kβkeCon,k,
where CDG is the installed capacity of the DG plant; αi is the ratio between exergy and energy content for fuel i which can be collected from relevant studies; qi is the energy content of per unit of fuel i; aCon,i is the amount of fuel j consumed to construct per capacity of the DG plant; xj is the chemical exergy content of the non-fuel material j; bCon,j is the mass of the non-fuel material j needed for per capacity of the DG plant; eCon,k is the emission amounts of pollutant k in the construction stage and the operation stage, respectively; βk is the abatement exergy factor of pollutant k.For the operation phase, the input of non-fuel materials is negligible, and thus the cumulative exergy consumption is mainly caused by fuel consumption which is proportional to power generation, shown as (5). For the DG plants that solely use renewable energies, the CECOpr can be neglected. Similarly, the emission abatement exergy of fuel-driven DG can be calculated by multiplying the amounts of different emissions and abatement exergy factors as shown in (6) and that of the renewable energy-driven DG is neglected,
(5)CECOpr=CDGHDGTDG∑iαiqiaOpr,iηOpr,i,(6)AbatEOpr=CDGHDGTDG∑kβkeOpr,k,
where HDG and TDG are the annual utilization hours (AUH) and the operation expectancy of the DG plant, respectively; aOpr,i is the cumulative mass of fuel j consumed to generate per unit of electric power; ηOpr,i is the overall efficiency of production and delivery for fuel i. eOpr,k is the emission amount of pollutant k in the construction stage and the operation stage, respectively.The cumulative exergy consumption and emission abatement exergy in the disposal/recycle phase consist of exergy consumption for unit disposal and exergy saving for unit recycle. In general, these exergies can be calculated proportionally to the exergies in the construction phase, shown as follows:(7)CECDis=-λCECCECCon,AbatEDis=-λAbatAbatECon,
where λCEC and λAbat represent the exergy saving ratios with respect to CECCon and AbatECon, respectively.


## 2.2. Optimization Model

The optimization variables of the siting and sizing problem are the locations of the DG plants in the studied distribution network and their installed capacities. The location variables can be expressed with integer numbers. Usually a DG unit is produced with a standard capacity, and thus the capacity variables of a DG plant can be expressed as integers that represent the number of installed DG units. Based on such arguments, the siting and sizing problem is modeled as a multi-objective integer programming problem in this paper. The objectives include minimizing the economic cost, life-cycle exergy consumption, and power losses.


### 2.2.1. Objective Functions

The multi-objective function can be written as(8)MinF={F1,F2,F3},
where F1 represents the life-cycle economic cost of utilizing DG plants; F2 represents the life-cycle exergy of using DG plants; F3 represents the expected power loss considering the uncertainties of power loads and DG outputs.The life-cycle economic cost of utilizing DG plants can be calculated as(9)F1=∑i∑j(CCon,ij+COpr,ij+CDis,ij),
where CCon,ij, COpr,ij and CDis,ij, are economic costs spent on construction, operation, and disposal/recycle to use jth type of DG at node i which can be calculated as
(10)CCon,ij=cCon,jCDG,ij,COpr,ij=ξCPV,jcOpr,jCDG,ijHDG,j,CDis,i=-ζTDG,j+1dCost,jCCon,ij,
where cCon,j is the economic cost for constructing per capacity of jth-type DG; cOpr,j is the economic cost for generating per unit of electricity by jth-type DG; dCost,j is the recycle coefficient for jth-type DG; coefficients ξCPV,j and ζ can be calculated as (11), in which f and δ stand for the inflation rate and discount rate, respectively. Generally, due to the intermittency of renewable energies, both the operation costs and power generation of DG plants that solely use renewable energies are much less than those of the DG plants which have the same capacities but use fuels:
(11)ξCPV,j=ζ-ζTDG,j+11-ζ,ζ=1+f1+δ.The objective function of minimizing life-cycle exergy can be written as(12)F2=∑i∑j(LCEDG,ij-ENoDG,ij),
where LCEDG,ij is the life-cycle exergy used by the DG plant of jth type at node i; ENoDG,ij stands for the expected exergy consumption to provide the same amount of DG-generated energies (electricity, heat, refrigeration, etc.) by alternative sources if DG plant ij is not installed.E
N
o
DG
,
i
j can be calculated as
(13)ENoDG,ij=CDG,ijHDG,j∑t=1TDG,j(γE,t+∑kυjkγk,t),
where γE,t represents the exergy needed for providing per unit of electric power in tth year by alternative sources; υjk is the amount of by-product k when generating per unit of electric power; γk,t stands for the exergy needed for providing per unit of by-product k in tth year by alternative sources.The objective function of minimizing expected power loss as(14)F3=E{ELoss},
where ELoss stands for the power loss and E{ELoss} represents the mean value of ELoss.


### 2.2.2. Constraints

The total capacity of installed DG plants should be no higher than a given cap for the sake of safety and reliability of power supply, that is,(15)∑i∑jCDG,ij≤CDGmax.Other main constraints include the nodal voltage magnitude limitations and the branch capacity limitations. In real operation, these limitations can be violated for short terms under some operating conditions. Considering the uncertainties of renewable energy sources and power loads, chance constraints are generally used to guarantee the practicability of the solution shown as follows:(16)Pr{Vimin<Vi<Vimax}≥ηV,Pr{|Si,j|<Si,jmax}≥ηS,
where Si,j stands for the apparent power flow in the branch from node i to node k; Vimin, Vimax, and Si,jmax are limitation values; ηV and ηS are required confidence parameters.For each operation condition, the load flow equations need to be obeyed, that is,(17)Pi=Vi∑j∈iVj[Gijcos(θi-θj)+Bijsin(θi-θj)],Qi=Vi∑j∈iVj[Gijsin(θi-θj)-Bijcos(θi-θj)],
where j∈i represents the set of every node j which is directly connected with node i by a branch in the distribution network, including node i itself; Pi and Qi stand for the active and reactive powers injected into node i; Vi and θi stand for the voltage magnitude and angle of node i; Gij and Bij represent the real and imaginary parts of admittance between node i and j.


## 2.2.1. Objective Functions

The multi-objective function can be written as(8)MinF={F1,F2,F3},
where F1 represents the life-cycle economic cost of utilizing DG plants; F2 represents the life-cycle exergy of using DG plants; F3 represents the expected power loss considering the uncertainties of power loads and DG outputs.The life-cycle economic cost of utilizing DG plants can be calculated as(9)F1=∑i∑j(CCon,ij+COpr,ij+CDis,ij),
where CCon,ij, COpr,ij and CDis,ij, are economic costs spent on construction, operation, and disposal/recycle to use jth type of DG at node i which can be calculated as
(10)CCon,ij=cCon,jCDG,ij,COpr,ij=ξCPV,jcOpr,jCDG,ijHDG,j,CDis,i=-ζTDG,j+1dCost,jCCon,ij,
where cCon,j is the economic cost for constructing per capacity of jth-type DG; cOpr,j is the economic cost for generating per unit of electricity by jth-type DG; dCost,j is the recycle coefficient for jth-type DG; coefficients ξCPV,j and ζ can be calculated as (11), in which f and δ stand for the inflation rate and discount rate, respectively. Generally, due to the intermittency of renewable energies, both the operation costs and power generation of DG plants that solely use renewable energies are much less than those of the DG plants which have the same capacities but use fuels:
(11)ξCPV,j=ζ-ζTDG,j+11-ζ,ζ=1+f1+δ.The objective function of minimizing life-cycle exergy can be written as(12)F2=∑i∑j(LCEDG,ij-ENoDG,ij),
where LCEDG,ij is the life-cycle exergy used by the DG plant of jth type at node i; ENoDG,ij stands for the expected exergy consumption to provide the same amount of DG-generated energies (electricity, heat, refrigeration, etc.) by alternative sources if DG plant ij is not installed.E
N
o
DG
,
i
j can be calculated as
(13)ENoDG,ij=CDG,ijHDG,j∑t=1TDG,j(γE,t+∑kυjkγk,t),
where γE,t represents the exergy needed for providing per unit of electric power in tth year by alternative sources; υjk is the amount of by-product k when generating per unit of electric power; γk,t stands for the exergy needed for providing per unit of by-product k in tth year by alternative sources.The objective function of minimizing expected power loss as(14)F3=E{ELoss},
where ELoss stands for the power loss and E{ELoss} represents the mean value of ELoss.


## 2.2.2. Constraints

The total capacity of installed DG plants should be no higher than a given cap for the sake of safety and reliability of power supply, that is,(15)∑i∑jCDG,ij≤CDGmax.Other main constraints include the nodal voltage magnitude limitations and the branch capacity limitations. In real operation, these limitations can be violated for short terms under some operating conditions. Considering the uncertainties of renewable energy sources and power loads, chance constraints are generally used to guarantee the practicability of the solution shown as follows:(16)Pr{Vimin<Vi<Vimax}≥ηV,Pr{|Si,j|<Si,jmax}≥ηS,
where Si,j stands for the apparent power flow in the branch from node i to node k; Vimin, Vimax, and Si,jmax are limitation values; ηV and ηS are required confidence parameters.For each operation condition, the load flow equations need to be obeyed, that is,(17)Pi=Vi∑j∈iVj[Gijcos(θi-θj)+Bijsin(θi-θj)],Qi=Vi∑j∈iVj[Gijsin(θi-θj)-Bijcos(θi-θj)],
where j∈i represents the set of every node j which is directly connected with node i by a branch in the distribution network, including node i itself; Pi and Qi stand for the active and reactive powers injected into node i; Vi and θi stand for the voltage magnitude and angle of node i; Gij and Bij represent the real and imaginary parts of admittance between node i and j.


## 2.3. Treatment of Uncertainties

### 2.3.1. Models of Uncertain Variables

Wind Turbine Output. Generally, the stochastic behaviors of the wind speed v in a wind farm can be described by the Weibull distribution whose probability density function (PDF) can be written as [37]
(18)fv(v)=kckvk-1exp[-(vc)k],
where k and c are shape factor and scale factor of the wind speed PDF at the wind farm.Ignoring minor nonlinearities, a turbine’s power output can be calculated with a given wind speedv:
(19)Pw={0v<vIorv>vOPwRv-vIvR-vIvI≤v≤vRPwRvR<v≤vO,
where vI, vO and vR are cut-in, cut-out, and rated wind speeds for the wind turbine; PwR is the turbine’s rated output.PV Output. Ignoring the minor factors, the output of a PV unit can be characterized as a function of solar irradiance and PV cell temperature [29, 36]:
(20)Ps=rrmax(1-κ·ΔT)Ps,max,
where r and rmax stand for the solar irradiance and its maximum value; Ps,max represents the maximum output at the rated cell temperature; ΔT is the deviation of the PV cell’s temperature from the rated one; κ is a coefficient.In the daytime, the PDF of solar irradiance can be characterized by Beta distribution [18, 29] which can be written as
(21)fr(r)=Γ(αr+βr)Γ(αr)Γ(βr)(rrmax)αr-1(1-rrmax)βr-1,
where αr and βr are distribution coefficients.The PV cell temperature is determined by the ambient temperature and prevailing wind speed. The deviation of the PV cell’s temperature can be described by normal distribution [29], that is,
(22)ΔT~N(0,σΔT2).Power Load. Generally, the uncertainty of the power load at a node can be described with normal distribution while the power factor can be regarded as a constant for simplicity [7, 38]. Therefore,
(23)PL,i~N(μL,i,σL,i2),
where PL,i represents the active power load at node i; μL,i and σL,i2 are the mean value and variance of PL,i.LHS-Based Semilinearized Probabilistic Load Flow Analysis. The typical procedure of LHS-based probabilistic load flow analysis can be found in [31], which consists of 3 steps: sampling, simulating, and accounting. In the sampling step, a required number of sample arrays of the uncertain variables are generated according to the variables’ probability density functions (see Section 2.3.1). Correlations of variables can be taken into account by Nataf transformation [31]. Each sample array represents a possible operating point and is used to evaluate the corresponding node voltages and power flows by solving (17) in the simulating step. Accounting on simulation results is then performed to obtain the statistical properties of certain variables, such as the probabilities in (16) and the expected power loss in (14).For simplicity, the load flow equation, that is, (17), can be expressed by (24) in which vector Y expresses the injected powers and vector X stands for the nodal voltages. Then, the power flows in branches can be expressed as a function of X.  Given an initial operating point Y0, the corresponding X0 and S0 can be calculated. However, since time-consuming iterations are needed due to the complexity and nonlinearity of (24), that is, (17), calculating X0 needs much more time than calculating S0 by (25), and this makes the simulating step requires much more time than the other 2 steps in PLF analysis:
(24)Y=g(X),(25)S=h(X).Expanding (24) at Y0 according to Taylor series and omitting the terms higher than the second order will yield
(26)X=X0+ΔX=X0+M0ΔY,
where vector ΔY expresses the changes of Y from Y0; vector ΔX stands for the changes of X; M0=J0-1, and J0 is the final Jacobian matrix at the initial operating point Y0. Equation (26) is also called “linearized load flow equation.”Similarly, (25) can be linearized as
(27)S=S0+ΔS=S0+T0ΔY,
where T0=L0M0 and L0=(∂h(X))/∂X|X=X0.Apparently, load flow simulation based on (26) and (27) will need much less calculating time than based on (24) and (25), and this is the theoretical fundament of cumulant-based series expansion methods for PLF analysis [29]. However, simulations show that (26) generates quite accurate results even with large values of ΔY; while, on the contrary, (27) can not always produce accurate results since it needs two times linearization from ΔY to ΔX. Nevertheless, (26) can be utilized to improve the calculation efficiency of the simulating step in LHS-based PLF analysis following a semi-linearized simulation method as below.Step 1.
Select the operating condition when all the uncertain variables are fixed at their mean values as the initial operating pointY0. Substitute Y0 into (24) and solve the corresponding load flow equations to get the initial nodal voltage vector X0 and the Jacobian matrix J0.Step 2.
Seti=1, and calculate ΔYi=Yi-Y0, where Yi is the ith sample array generated by LHS method in the sampling step.Step 3.
SubstituteΔYi into (26) and calculate Xi.Step 4.
SubstituteXi into (25) to calculate Si.Step 5.
Seti=i+1, repeat Step 2 to Step 5 until all the sample arrays are evaluated.Thus, the semi-linearized simulation method avoids the time-consuming iterations in calculating the nodal voltages for different operating condition samples and guarantees the accuracy by keeping the nonlinearity of (25). Meanwhile, by using LHS-based PLF analysis framework, the inaccuracy and instability of cumulant-based series expansion are avoided in the PLF analysis.


## 2.3.1. Models of Uncertain Variables

Wind Turbine Output. Generally, the stochastic behaviors of the wind speed v in a wind farm can be described by the Weibull distribution whose probability density function (PDF) can be written as [37]
(18)fv(v)=kckvk-1exp[-(vc)k],
where k and c are shape factor and scale factor of the wind speed PDF at the wind farm.Ignoring minor nonlinearities, a turbine’s power output can be calculated with a given wind speedv:
(19)Pw={0v<vIorv>vOPwRv-vIvR-vIvI≤v≤vRPwRvR<v≤vO,
where vI, vO and vR are cut-in, cut-out, and rated wind speeds for the wind turbine; PwR is the turbine’s rated output.PV Output. Ignoring the minor factors, the output of a PV unit can be characterized as a function of solar irradiance and PV cell temperature [29, 36]:
(20)Ps=rrmax(1-κ·ΔT)Ps,max,
where r and rmax stand for the solar irradiance and its maximum value; Ps,max represents the maximum output at the rated cell temperature; ΔT is the deviation of the PV cell’s temperature from the rated one; κ is a coefficient.In the daytime, the PDF of solar irradiance can be characterized by Beta distribution [18, 29] which can be written as
(21)fr(r)=Γ(αr+βr)Γ(αr)Γ(βr)(rrmax)αr-1(1-rrmax)βr-1,
where αr and βr are distribution coefficients.The PV cell temperature is determined by the ambient temperature and prevailing wind speed. The deviation of the PV cell’s temperature can be described by normal distribution [29], that is,
(22)ΔT~N(0,σΔT2).Power Load. Generally, the uncertainty of the power load at a node can be described with normal distribution while the power factor can be regarded as a constant for simplicity [7, 38]. Therefore,
(23)PL,i~N(μL,i,σL,i2),
where PL,i represents the active power load at node i; μL,i and σL,i2 are the mean value and variance of PL,i.LHS-Based Semilinearized Probabilistic Load Flow Analysis. The typical procedure of LHS-based probabilistic load flow analysis can be found in [31], which consists of 3 steps: sampling, simulating, and accounting. In the sampling step, a required number of sample arrays of the uncertain variables are generated according to the variables’ probability density functions (see Section 2.3.1). Correlations of variables can be taken into account by Nataf transformation [31]. Each sample array represents a possible operating point and is used to evaluate the corresponding node voltages and power flows by solving (17) in the simulating step. Accounting on simulation results is then performed to obtain the statistical properties of certain variables, such as the probabilities in (16) and the expected power loss in (14).For simplicity, the load flow equation, that is, (17), can be expressed by (24) in which vector Y expresses the injected powers and vector X stands for the nodal voltages. Then, the power flows in branches can be expressed as a function of X.  Given an initial operating point Y0, the corresponding X0 and S0 can be calculated. However, since time-consuming iterations are needed due to the complexity and nonlinearity of (24), that is, (17), calculating X0 needs much more time than calculating S0 by (25), and this makes the simulating step requires much more time than the other 2 steps in PLF analysis:
(24)Y=g(X),(25)S=h(X).Expanding (24) at Y0 according to Taylor series and omitting the terms higher than the second order will yield
(26)X=X0+ΔX=X0+M0ΔY,
where vector ΔY expresses the changes of Y from Y0; vector ΔX stands for the changes of X; M0=J0-1, and J0 is the final Jacobian matrix at the initial operating point Y0. Equation (26) is also called “linearized load flow equation.”Similarly, (25) can be linearized as
(27)S=S0+ΔS=S0+T0ΔY,
where T0=L0M0 and L0=(∂h(X))/∂X|X=X0.Apparently, load flow simulation based on (26) and (27) will need much less calculating time than based on (24) and (25), and this is the theoretical fundament of cumulant-based series expansion methods for PLF analysis [29]. However, simulations show that (26) generates quite accurate results even with large values of ΔY; while, on the contrary, (27) can not always produce accurate results since it needs two times linearization from ΔY to ΔX. Nevertheless, (26) can be utilized to improve the calculation efficiency of the simulating step in LHS-based PLF analysis following a semi-linearized simulation method as below.Step 1.
Select the operating condition when all the uncertain variables are fixed at their mean values as the initial operating pointY0. Substitute Y0 into (24) and solve the corresponding load flow equations to get the initial nodal voltage vector X0 and the Jacobian matrix J0.Step 2.
Seti=1, and calculate ΔYi=Yi-Y0, where Yi is the ith sample array generated by LHS method in the sampling step.Step 3.
SubstituteΔYi into (26) and calculate Xi.Step 4.
SubstituteXi into (25) to calculate Si.Step 5.
Seti=i+1, repeat Step 2 to Step 5 until all the sample arrays are evaluated.Thus, the semi-linearized simulation method avoids the time-consuming iterations in calculating the nodal voltages for different operating condition samples and guarantees the accuracy by keeping the nonlinearity of (25). Meanwhile, by using LHS-based PLF analysis framework, the inaccuracy and instability of cumulant-based series expansion are avoided in the PLF analysis.


## 2.4. Solving Approach and Methodology Framework

To solve the proposed multi-objective optimization model, the nondominated sorting genetic algorithm II (NSGA-II) [39] is used in this paper. NSGA-II greatly improves the performances of NSGA by introducing elitism and diversity-preserving mechanism algorithms and assures better convergence without losing solution diversity [40]. It should be noted that, instead of giving one optimal solution, NSGA and NSGA-II generate a set of nondominated solutions, that is, Pareto-optimal front, for multi-objective optimization problems according to the Pareto optimality principle [41]. Each solution actually represents the optimal solution under a different tradeoff relationship of the objectives. This enables the planner to see different possible solutions and helps the planner to choose a final solution based on certain tradeoff criteria according to a specific situation after comparing possible alternatives. The brief framework of the optimization algorithm used in this paper is shown in Figure 1, where the step of evaluating chromosomes consists of several substeps shown in Figure 2.Figure 1
Brief framework of the NSGA-II algorithm in this paper.Figure 2
Substeps in the step of evaluating chromosomes.Each chromosome in the algorithm represents a candidate solution for the optimization problem and is coded with integer genes as shown in Figure3, where each pane represents a gene in the chromosome. The position of the gene indicates the candidate site and type of a DG plant, and the value of the gene stands for the number of installed units.Figure 3
Coding scheme of the genes in a chromosome.


## 3. Case Study

### 3.1. System Data

The system under study is the typical IEEE (Institute of Electrical and Electronics Engineers) 33-node system as shown in Figure4. The network parameters are available from [42, 43]. Candidate DG units include wind turbines (WT), photovoltaic units (PV), and micronatural gas turbines (MNGT). Possible sites and sizes of DG plants are listed in Table 1, in which each integer number of a DG plant represent the feasible maximum number of installed DG units. Each DG unit is assumed to have a rated power output of 20 kW. DG units of the same type have the same technical parameters regardless of the different sites. The cut-in, rated, and cut-out wind speeds for wind turbines are 4 m/s, 14 m/s, and 20 m/s, respectively. MNGT units are used as combined heat and power generators, which produce 1 kWh useful heat as a by-product when generating 1 kWh electricity. For simplicity, the statistical characteristics of wind speeds, solar irradiances, and PV cell temperatures are assumed to be, respectively, the same for different sites considering the short distances between them, which are listed in Table 2. For the power loads, the mean values are the constant values in [42], and the standard deviations are set as 5% of the corresponding mean values. Randomly generated positive definite matrices are used as the correlation matrices of uncertain variables [31]. The exergy consumption and economic cost parameter are collected from [27, 44] and listed in Table 3. The exergy cost for supplying 1 kWh electricity and 1 kWh heat by alternative sources are 19, 162 kJ and 10, 123 kJ, respectively, regardless of possible changes in the future. Both the required confidences for nodal voltage and branch apparent power flow (ηV and ηS) are 95%. The inflation rate and discount rate are 3% and 5.35%, respectively.Table 1
Candidate sites and sizes of DG plants.
Node
16
17
18
21
22
23
24
25
31
32
Total
WT
4
4
4
3
3
3
4
3
4
4
36
PV
4
4
4
3
4
3
4
3
3
3
35
MNGT
2
2
2
0
3
2
2
3
2
0
18Table 2
Parameters of wind speed, solar irradiance, and PV cell temperature [7, 29, 36].
k
c
α
r
β
r
κ
σ
Δ
T
2
6.5
5.5
1.8
0.0035
5Table 3
Cost and operation parameters of DG units.
c
Con ($/kW)
c
Opr ($/kWh)
d
Cost
T
DG (yr)
H
DG (h)
CEC (GJ)
AbatE (GJ)
E
N
o
DG (GJ)
WT
1,800
0.05
0.05
25
2,100
118.3
50.8
20,120.1
PV
2,000
0.03
0.05
30
1,650
406.5
189.2
18,970.4
MNGT
850
0.04
0.05
30
6,000
33,348.7
9,701.2
105,426.7Figure 4
Structure of IEEE 33-node test system and candidate DG sites.The optimization problem contains 28 integer variables, and the maximum number of uncertain variables in PLF reaches 52. The population of chromosomes in NSGA-II is set as 200 and the maximum generation of evolution as 100, which are large enough to guarantee the convergence of the algorithm.


### 3.2. Results and Discussions

To verify the proposed semi-linearized simulation method, probabilistic load flow analyses of 100 randomly-generated solutions are carried out first. Figures5(a) and 5(b) show the comparison of the results calculated by the proposed method and the original LHS-based PLF method, where Pr{Si,j} represents the overload probability of the branch from node i to nodej, that is, Pr{Si,j} = 1-Pr{|Si,j|<Si,jmax}. For simplicity, only the results of the branch from node 14 to node 15 and the branch from node 30 to node 31 are pictured since results of other branches are quite similar. It can be seen that the proposed method generates results very close to the original method. The largest difference is about 0.0035 which can be neglected in analyses. Figure 5(c) shows the difference in calculation time of the two methods. In average, the proposed method uses about 0.072 s to perform a PLF analysis while the original method needs 0.941 s. That is to say, the proposed method can save about 92.3% of the calculation time for PLF analysis. These results demonstrate the accuracy and efficiency of the proposed LHS-based semi-linearized PLF analysis method.Comparison of the results from the proposed method and the original method.
(a)
(b)
(c)Figure6 shows the final Pareto-optimal front of the proposed optimization model. The optimization results indicate that as the economic cost of DG plants increases, that is, as the installed capacity of DG plants grows, the cost of exergy and expected power loss become lower. Since it costs more exergy to supply electricity and heat by other alternative sources than by DG plants, the exergy cost of installing DG plants is negative. The absolute value of exergy represents the saved exergy and indicates the benefit for maintaining the sustainability of using the DG plants in this case. If entirely converted by the equivalent exergy contents, 1 PJ of the saved exergy means a reduction of 45, 126.4 tons of coal consumption or 170, 648.5 tons of CO2 emission. The maximum exergy that can be saved reaches 2.33 PJ, and the economic cost for this is 10.33 million US dollars. Before integrating DG plants, the expected power loss of the system reaches 202.9 kW. As the installed capacity of DG increases, the expected power loss can be reduced down to 133.5 kW, saving about 34.2%. Thus, it can be concluded that integrating DG plants properly into the distribution system can efficiently improve the system’s sustainability and reduce the power loss.Pareto-optimal front of the optimization result. (a) The 3-dimensional distribution of the Pareto-optimal front. (b) The cost-exergy distribution of the Pareto-optimal front. (c) The cost-loss distribution of the Pareto-optimal front.
(a)
(b)
(c)Figures6(b) and 6(c) show the changes in tradeoff relationships of the 3 objective functions in the Pareto-optimal front. Since the exergy saving of installing a certain DG unit is constant, the tradeoff relationship between the economic cost and saved exergy is close to linear relation. That is to say, if the planner wants to save a certain amount of exergy, a certain amount of economic cost is needed. Comparatively, the tradeoff relationship between the economic cost and the power loss is nonlinear. As the economic cost increases, the reduction in power loss caused by adding per unit of economic cost gets smaller. This indicates that adding more DG plants to reduce power loss is not so effective when there are already plenty of them.The Pareto-optimal front of the proposed optimization model provides useful information for the system planner to set a desirable arrangement for the DG plants. It is convenient to choose one of the Pareto-optimal solutions as the final decision if the planner is satisfied with the corresponding tradeoff relationship according to a specific situation. Although usually not all of the Pareto-optimal solutions are provided in optimization, the distribution of these solutions helps the planner to investigate the possible range of the costs and benefits of the solutions, and thus the planner can determine the optimal tradeoff among the different objectives and at the same time set rational expectations for them. For example, if investment budget is limited, the planner may prefer solutions with lower economic costs but at the same time he or she has to accept the smaller benefits for improving environmental sustainability and reducing power loss.Figure7 shows the contributions and the variation tendencies of the three types of DG plants in the total DG capacity, economic cost, and saved exergy as the total economic cost increases in the Pareto-optimal front. All the installed capacities of WT, PV, and MNGT tend to increase when the total economic cost increases. The installed capacity of MNGT is always smaller than those of WT and PV but has a larger increase rate. The average installed capacities of WT, PV, and MNGT are 309.0 kW, 546.7 kW and 20 kW respectively when the total economic cost is lower than 5 million US dollars, which reach 644.7 kW, 675.0 kW and 250.9 kW, respectively, when the total economic cost is higher than 7 million US dollars. This indicates that MNGT is not the prior choice compared with WT, and PV when the economic cost is limited. Although an MNGT unit can bring 2.1~2.3 fold more exergy saving than WT unit, and PV unit, the much higher economic cost, which is 3.7~4.9 fold higher, prevents its application. The same amount of exergy saving caused by installing MNGT a unit can be obtained by installing more WT or PV units with lower cost if the constraints can be satisfied. However, as a reliable and stable power source, the MNGT units can be used to suppress the negative impacts of WT and PV units on the distribution system, and this is why the MNGT units can not be replaced by other DG units in the solutions of Pareto-optimal front.Contributions of WT, PV, and MNGT to the total DG capacity, economic cost, and saved exergy.
(a)
(b)
(c)Since the system is already close to its upper limit of power distribution due to the small capacities of branches, the uncertainties of power loads bring great overload risks. The assumption of tight correlations between the power loads at different nodes further increases the overload risks as the power loads tend to vary unanimously. Figure6 shows the overload probabilities of distribution branches before and after the DG plants are installed. The applied siting and sizing plan of DG plants in Figure 6 is the minimum economic cost plan, whose details are shown in Table 4. This plan costs 1.60 million US dollars, −0.47 PJ of exergy and 174.1 kW of expected power loss. The overload probabilities of branches are obviously reduced by integrating DG plants as shown in Figure 8. Additionally, simulations show that more expensive plans generally perform better in reducing the overload risks than this plan.Table 4
Installed units of the minimum economic cost solution.
Node
16
17
18
21
22
23
24
25
31
32
WT
1
1
0
0
0
1
2
1
1
1
PV
4
3
2
0
1
0
2
2
0
3
MNGT
0
0
0
0
0
0
0
0
0
0Figure 8
Overload probabilities of distribution branches. The dash line represents the overload probability of 5%.


## 3.1. System Data

The system under study is the typical IEEE (Institute of Electrical and Electronics Engineers) 33-node system as shown in Figure4. The network parameters are available from [42, 43]. Candidate DG units include wind turbines (WT), photovoltaic units (PV), and micronatural gas turbines (MNGT). Possible sites and sizes of DG plants are listed in Table 1, in which each integer number of a DG plant represent the feasible maximum number of installed DG units. Each DG unit is assumed to have a rated power output of 20 kW. DG units of the same type have the same technical parameters regardless of the different sites. The cut-in, rated, and cut-out wind speeds for wind turbines are 4 m/s, 14 m/s, and 20 m/s, respectively. MNGT units are used as combined heat and power generators, which produce 1 kWh useful heat as a by-product when generating 1 kWh electricity. For simplicity, the statistical characteristics of wind speeds, solar irradiances, and PV cell temperatures are assumed to be, respectively, the same for different sites considering the short distances between them, which are listed in Table 2. For the power loads, the mean values are the constant values in [42], and the standard deviations are set as 5% of the corresponding mean values. Randomly generated positive definite matrices are used as the correlation matrices of uncertain variables [31]. The exergy consumption and economic cost parameter are collected from [27, 44] and listed in Table 3. The exergy cost for supplying 1 kWh electricity and 1 kWh heat by alternative sources are 19, 162 kJ and 10, 123 kJ, respectively, regardless of possible changes in the future. Both the required confidences for nodal voltage and branch apparent power flow (ηV and ηS) are 95%. The inflation rate and discount rate are 3% and 5.35%, respectively.Table 1
Candidate sites and sizes of DG plants.
Node
16
17
18
21
22
23
24
25
31
32
Total
WT
4
4
4
3
3
3
4
3
4
4
36
PV
4
4
4
3
4
3
4
3
3
3
35
MNGT
2
2
2
0
3
2
2
3
2
0
18Table 2
Parameters of wind speed, solar irradiance, and PV cell temperature [7, 29, 36].
k
c
α
r
β
r
κ
σ
Δ
T
2
6.5
5.5
1.8
0.0035
5Table 3
Cost and operation parameters of DG units.
c
Con ($/kW)
c
Opr ($/kWh)
d
Cost
T
DG (yr)
H
DG (h)
CEC (GJ)
AbatE (GJ)
E
N
o
DG (GJ)
WT
1,800
0.05
0.05
25
2,100
118.3
50.8
20,120.1
PV
2,000
0.03
0.05
30
1,650
406.5
189.2
18,970.4
MNGT
850
0.04
0.05
30
6,000
33,348.7
9,701.2
105,426.7Figure 4
Structure of IEEE 33-node test system and candidate DG sites.The optimization problem contains 28 integer variables, and the maximum number of uncertain variables in PLF reaches 52. The population of chromosomes in NSGA-II is set as 200 and the maximum generation of evolution as 100, which are large enough to guarantee the convergence of the algorithm.


## 3.2. Results and Discussions

To verify the proposed semi-linearized simulation method, probabilistic load flow analyses of 100 randomly-generated solutions are carried out first. Figures5(a) and 5(b) show the comparison of the results calculated by the proposed method and the original LHS-based PLF method, where Pr{Si,j} represents the overload probability of the branch from node i to nodej, that is, Pr{Si,j} = 1-Pr{|Si,j|<Si,jmax}. For simplicity, only the results of the branch from node 14 to node 15 and the branch from node 30 to node 31 are pictured since results of other branches are quite similar. It can be seen that the proposed method generates results very close to the original method. The largest difference is about 0.0035 which can be neglected in analyses. Figure 5(c) shows the difference in calculation time of the two methods. In average, the proposed method uses about 0.072 s to perform a PLF analysis while the original method needs 0.941 s. That is to say, the proposed method can save about 92.3% of the calculation time for PLF analysis. These results demonstrate the accuracy and efficiency of the proposed LHS-based semi-linearized PLF analysis method.Comparison of the results from the proposed method and the original method.
(a)
(b)
(c)Figure6 shows the final Pareto-optimal front of the proposed optimization model. The optimization results indicate that as the economic cost of DG plants increases, that is, as the installed capacity of DG plants grows, the cost of exergy and expected power loss become lower. Since it costs more exergy to supply electricity and heat by other alternative sources than by DG plants, the exergy cost of installing DG plants is negative. The absolute value of exergy represents the saved exergy and indicates the benefit for maintaining the sustainability of using the DG plants in this case. If entirely converted by the equivalent exergy contents, 1 PJ of the saved exergy means a reduction of 45, 126.4 tons of coal consumption or 170, 648.5 tons of CO2 emission. The maximum exergy that can be saved reaches 2.33 PJ, and the economic cost for this is 10.33 million US dollars. Before integrating DG plants, the expected power loss of the system reaches 202.9 kW. As the installed capacity of DG increases, the expected power loss can be reduced down to 133.5 kW, saving about 34.2%. Thus, it can be concluded that integrating DG plants properly into the distribution system can efficiently improve the system’s sustainability and reduce the power loss.Pareto-optimal front of the optimization result. (a) The 3-dimensional distribution of the Pareto-optimal front. (b) The cost-exergy distribution of the Pareto-optimal front. (c) The cost-loss distribution of the Pareto-optimal front.
(a)
(b)
(c)Figures6(b) and 6(c) show the changes in tradeoff relationships of the 3 objective functions in the Pareto-optimal front. Since the exergy saving of installing a certain DG unit is constant, the tradeoff relationship between the economic cost and saved exergy is close to linear relation. That is to say, if the planner wants to save a certain amount of exergy, a certain amount of economic cost is needed. Comparatively, the tradeoff relationship between the economic cost and the power loss is nonlinear. As the economic cost increases, the reduction in power loss caused by adding per unit of economic cost gets smaller. This indicates that adding more DG plants to reduce power loss is not so effective when there are already plenty of them.The Pareto-optimal front of the proposed optimization model provides useful information for the system planner to set a desirable arrangement for the DG plants. It is convenient to choose one of the Pareto-optimal solutions as the final decision if the planner is satisfied with the corresponding tradeoff relationship according to a specific situation. Although usually not all of the Pareto-optimal solutions are provided in optimization, the distribution of these solutions helps the planner to investigate the possible range of the costs and benefits of the solutions, and thus the planner can determine the optimal tradeoff among the different objectives and at the same time set rational expectations for them. For example, if investment budget is limited, the planner may prefer solutions with lower economic costs but at the same time he or she has to accept the smaller benefits for improving environmental sustainability and reducing power loss.Figure7 shows the contributions and the variation tendencies of the three types of DG plants in the total DG capacity, economic cost, and saved exergy as the total economic cost increases in the Pareto-optimal front. All the installed capacities of WT, PV, and MNGT tend to increase when the total economic cost increases. The installed capacity of MNGT is always smaller than those of WT and PV but has a larger increase rate. The average installed capacities of WT, PV, and MNGT are 309.0 kW, 546.7 kW and 20 kW respectively when the total economic cost is lower than 5 million US dollars, which reach 644.7 kW, 675.0 kW and 250.9 kW, respectively, when the total economic cost is higher than 7 million US dollars. This indicates that MNGT is not the prior choice compared with WT, and PV when the economic cost is limited. Although an MNGT unit can bring 2.1~2.3 fold more exergy saving than WT unit, and PV unit, the much higher economic cost, which is 3.7~4.9 fold higher, prevents its application. The same amount of exergy saving caused by installing MNGT a unit can be obtained by installing more WT or PV units with lower cost if the constraints can be satisfied. However, as a reliable and stable power source, the MNGT units can be used to suppress the negative impacts of WT and PV units on the distribution system, and this is why the MNGT units can not be replaced by other DG units in the solutions of Pareto-optimal front.Contributions of WT, PV, and MNGT to the total DG capacity, economic cost, and saved exergy.
(a)
(b)
(c)Since the system is already close to its upper limit of power distribution due to the small capacities of branches, the uncertainties of power loads bring great overload risks. The assumption of tight correlations between the power loads at different nodes further increases the overload risks as the power loads tend to vary unanimously. Figure6 shows the overload probabilities of distribution branches before and after the DG plants are installed. The applied siting and sizing plan of DG plants in Figure 6 is the minimum economic cost plan, whose details are shown in Table 4. This plan costs 1.60 million US dollars, −0.47 PJ of exergy and 174.1 kW of expected power loss. The overload probabilities of branches are obviously reduced by integrating DG plants as shown in Figure 8. Additionally, simulations show that more expensive plans generally perform better in reducing the overload risks than this plan.Table 4
Installed units of the minimum economic cost solution.
Node
16
17
18
21
22
23
24
25
31
32
WT
1
1
0
0
0
1
2
1
1
1
PV
4
3
2
0
1
0
2
2
0
3
MNGT
0
0
0
0
0
0
0
0
0
0Figure 8
Overload probabilities of distribution branches. The dash line represents the overload probability of 5%.


## 4. Conclusion

In this paper, a sustainability-oriented multi-objective optimization model for siting and sizing DG plants in distribution systems is proposed. In the model, life cycle exergy is used as a unified indicator of all the impacts on environmental sustainability of integrating DG plants, and it is optimized in the model to obtain better environmental sustainability. Exergetic life cycle analysis methodology is used to evaluate the life cycle exergy of the planned system. To control the risks brought by the unstable DG outputs and power loads, chance constraints are used in the optimization model. Additionally, a semi-linearized simulation method is proposed to improve the efficiency of probabilistic load flow analysis based on Latin hypercube sampling.A case study based on modified IEEE 33-node system is carried out to demonstrate the proposed model and method. Numerical results show that the semi-linearized simulation method greatly improves the calculation efficiency of probabilistic load flow analysis, and the accuracy loss is little. The proposed optimization model generates a series of possible solutions known as the Pareto-optimal front and provides useful information for the system planner to set a desirable arrangement for the DG plants. The Pareto-optimal solutions in the case study indicate that with proper arrangement DG plants can produce considerable benefits for maintaining environmental sustainability, reducing power losses, and controlling operation risks. Wind turbines and photovoltaic units are prior choices compared to the micronatural gas turbines when the total economic cost is limited under a low level in this case study. More integration of micro natural gas turbines is beneficial for maintaining environmental sustainability and reducing power losses but obviously increases the economic cost.How to finally determine the sites and sizes of the DG plants according to the Pareto-optimal solutions is another important problem that needs methods for multicriteria decision making. Further study on this problem should be carried out in the future. Additionally, some modifications can be made based on the proposed model for studies such as optimization of reactive power compensation capacities and integration of electric vehicle charging stations.



---

*Source: 291930-2013-06-06.xml*