- Email: [email protected]
Optimal power flow solutions using differential evolution algorithm integrated with effective constraint handling techniques
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Contents lists available at ScienceDirect
Engineering Applications of Artificial Intelligence journal homepage: www.elsevier.com/locate/engappai
Optimal power flow solutions using differential evolution algorithm integrated with effective constraint handling techniques Partha P. Biswas a, *, P.N. Suganthan a , R. Mallipeddi b, *, Gehan A.J. Amaratunga c a b c
School of Electrical and Electronic Engineering, Nanyang Technological University, Singapore Kyungpook National University, Daegu, South Korea Department of Engineering, University of Cambridge, UK
a r t i c l e
i n f o
Keywords: Optimal power flow Power loss Emission Differential evolution (DE) Constraint handling (CH) Ensemble of constraint handling techniques (ECHT)
a b s t r a c t Optimal power flow (OPF) is a highly non-linear complex optimization problem where the steady state parameters of an electrical network need to be determined for its economical and efficient operation. The complexity of the problem escalates with ubiquitous presence of constraints in the problem. Solving OPF remains a popular but challenging task among power system researchers. In last couple of decades, numerous evolutionary algorithms (EAs) have been applied to find optimal solutions with different objectives of OPF. However, the search method adopted by EAs is unconstrained. An extensively used methodology to discard infeasible solutions found during the search process is the static penalty function approach. The process requires appropriate selection of penalty coefficients decided largely by tedious trial and error method. This paper presents performance evaluation of proper constraint handling (CH) techniques — superiority of feasibly solutions (SF), self-adaptive penalty (SP) and an ensemble of these two constraint handling techniques (ECHT) with differential evolution (DE) being the basic search algorithm, on the problem of OPF. The methods are tested on standard IEEE 30, IEEE 57 and IEEE 118-bus systems for several OPF objectives such as cost, emission, power loss, voltage stability etc. Single objective and weighted sum multi-objective cases of OPF are studied under the scope of this literature. Simulation results are analyzed and compared with most recent studies on the problem. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction Since its inception more than half a century ago, the optimal power flow (OPF) remains a widely cultivated topic among power system research communities across the globe due to the intriguing multifaceted challenges it poses. The OPF is formulated as single or multiobjective problem of minimizing fuel cost, emission, transmission loss, voltage deviation etc. with constraints to be satisfied on generator capability, line capacity, bus voltage and power flow balance. Output of OPF program helps to determine the optimal operating state of a power system and the corresponding settings of control variables for economic and secure operation (Cai et al., 2008). Major control variables refer to the generated real power and generator bus voltages of the network. The latter controls the reactive power flow which is further compensated by adding capacitor banks of appropriate ratings to the network feeding usually the inductive loads. Voltage vectors of load buses and complex power flows in the lines determined during the process of optimization
represent the system optimal operating state that would result in single or multiple objectives of the network being largely fulfilled. In summary, OPF involves intricate calculations with multiple variables and finding optimal solution while satisfying all constraints simultaneously is the most difficult part one encounters. In earlier days, in use were classical numerical optimization methods which suffered from convexity, assumption of continuity and normally employed a gradient based search that converged to local optima. Revolution in numerical optimization introduced several evolutionary algorithms and techniques in last few decades. Most of these methods can successfully overcome the problem of premature convergence and are able to explore the search space in pursuit of global optima. OPF problem has seen application of numerous such evolutionary algorithms. A few standard objectives of OPF were optimized in Abaci and Yamacli (2016) for IEEE bus systems using differential search algorithm (DSA), an effective algorithm for real-valued numerical optimization problems. Daryani et al. (2016) improved standard group
* Corresponding authors.
E-mail addresses: [email protected] (P.P. Biswas), [email protected] (P.N. Suganthan), [email protected] (R. Mallipeddi), [email protected] (G.A.J. Amaratunga). https://doi.org/10.1016/j.engappai.2017.10.019 Received 24 June 2017; Received in revised form 5 October 2017; Accepted 25 October 2017 0952-1976/© 2017 Elsevier Ltd. All rights reserved.
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
members would result in restricted scope of exploration and exploitation by the algorithm during the search process. System security constraints gained importance in determining maximum loading limits of selected networks in some literatures using genetic algorithm (GA) (Acharjee, 2012), fuzzy logic (Mallick et al., 2013) and chaotic PSO (Acharjee et al., 2011). In optimal reactive power flow problem, voltage security constraint is considered in Rabiee and Parniani (2013). However, state-of-the-art constraint handling (CH) techniques in OPF problem remain largely untested. In present study, superiority of feasible solutions (SF), self-adaptive penalty (SP) methods and an ensemble of these two constraint handling techniques (ECHT) are applied to OPF problem with differential evolution (DE) (Storn and Price, 1997) as the basic search algorithm. Performance of these CH techniques is assessed, results are statistically compared and analyzed in detail. It is worthwhile to note that OPF problem has both equality and inequality constraints. Equality constraints are the power balance equations where both active and reactive power generated in the network must be equal to the demand and losses in the network. Convergence of power flow to a solution ensures that the power balance equations are automatically satisfied. Inequality constraints on slack generator power, reactive power output of the generators, bus voltage limits and line capacities need special attention. SF and SP methods of constraint handling work differently on their own way in handling inequality constraints. The motivation behind ECHT is to assess performance of both the constraint handling techniques by the combinatorial algorithm (i.e. ECHT) as one technique might work better than the other on some problems but not on all problems. Instead of user finding the superior performer out of the two for a particular problem, ECHT takes the burden to apply both the techniques and output close to, if not the best results obtained by one of the constituting techniques. The constraint handling (CH) methods are all successfully implemented in OPF problem on standard IEEE 30, IEEE 57 and IEEE 118-bus systems in our research presented in this paper. Generation cost, emission, real power loss and voltage stability of the network are all individually (as single objective) optimized. Study cases with multi-fuel options for the generators and valve-point loading effect are also considered. Although the study cases mostly consider all continuous variables, handling of discrete variables is also discussed and a study-case results with discrete variables are appended to the result section. Weighted sum approach is adopted in multi-objective optimization for selected study cases of both 30-bus and 57-bus systems. Simulation results are compared and critically investigated against constraint violation with most recent studies on optimal power flow (OPF) found in the literatures. A summary of contributions of this work can be listed pointwise as below:
search optimization algorithm with adaptation to develop adaptive group search optimization (AGSO) to perform study on OPF. Chaib et al. (2016) applied backtracking search optimization algorithm (BSA) to perform OPF calculation with multi-fuel options and valve-point effect in thermal generators. Improved colliding bodies optimization (ICBO) algorithm (Bouchekara et al., 2016b) increased the number of colliding bodies in each iteration to enhance performance of the algorithm when applied to the problem of OPF. Mohamed et al. (2017) applied moth swarm algorithm (MSA) on numerous objectives of OPF for various bus systems to exhibit fast execution time and quick convergence of the algorithm. Chaos theory was incorporated in artificial bee colony to form chaotic artificial bee colony (CABC) in Ayan et al. (2015) and security constrained OPF was solved. Agent based gravitational search algorithm (GSA), where masses of agents play vital roles in guiding the search process, was adopted in Bhowmik and Chakraborty (2015) to solve OPF. To enhance exploration capability and population diversity of biogeography based optimization algorithm (BBO), adaptive real coded BBO (ARCBBO) was suggested in Kumar and Premalatha (2015) and applied to OPF problem. Dynamic adjustment of control parameters in adaptive partitioning flower pollination algorithm (APFPA) (Mahdad and Srairi, 2016) showed higher convergence speed and better accuracy than FPA for OPF solutions. In fuzzy harmony search algorithm (FHSA) (Pandiarajan and Babulal, 2016), effect of fuzzy based automatic adjustment of pitch adjustment rate and bandwidth of harmony search algorithm was studied in applying to the problem of OPF. Fuzzy based adaptive configuration of particle swarm optimization was established in Naderi et al. (2017) to solve a sub-problem of OPF. In solving OPF, Levy mutation strategy for teaching–learning based optimization (LTLBO) (Ghasemi et al., 2015) introduced a Levy mutation operator to enhance exploration at early search stages. Krill herd algorithm (KHA) (Roy and Paul, 2015) and its variant stud krill herd algorithm (SKH) (Pulluri et al., 2017, 2016) have also been popular in finding OPF solutions. Grenade explosion method (GEM) (Bouchekara et al., 2016a), glowworm swarm optimization (GSO) (Reddy and Rathnam, 2016) and Gbest guided artificial bee colony (GABC) (Roy and Jadhav, 2015) have also shown promising results for OPF as reported in respective literatures. Modified imperialist competitive algorithm (MICA) and teaching–learning algorithm (TLA) were hybridized (MICA-TLA) in Ghasemi et al. (2014) to improve local search and convergence of original ICA algorithm in finding OPF solutions. Multi-objective formulation and solution of OPF are found in Shaheen et al. (2016, 2017), Sedighizadeh et al. (2011) and Niknam et al. (2012). Our study focuses on single and weighted-sum multiobjective optimization cases of OPF. As evident from the literature review, either basic or the improved version of an algorithm was tried for real-parameter optimization problem of OPF. Each of the methods has its own strengths and weaknesses, and as depicted in No Free Lunch theorem (Wolpert and Macready, 1997), no single optimization method is capable of solving all types of real-world problems in an effective manner. Differential evolution (DE) and its variants are found to be among the top performing algorithms for single objective real parameter optimization as reflected in the results of CEC competitions (Liang et al., 2006; Li et al., 2008; Liang et al., 2013; Chen et al., 2014). We adopt DE as the basic search algorithm for optimization of OPF objectives. On the aspect constraint handling in OPF, literature (Daryani et al., 2016) defined security objective entwining a couple of inequality constraints as one of the multiple objectives to be minimized. Rest all literatures adopted either static penalty function approach or straightaway discarded the population members that led to constraint violation. As mentioned beforehand, penalty function approach is sensitive to selection of penalty coefficient. Small penalty coefficient over-explores the infeasible region, thus delaying the process of finding feasible solutions, and may prematurely converge to an infeasible solution. On the other hand, large penalty coefficient may not explore the infeasible region properly, thereby resulting in untimely convergence (Mallipeddi et al., 2012). In second approach, removal of infeasible population
∙ Applying CH techniques, superiority of feasible solutions (SF) and self-adaptive penalty (SP) individually to handle constraints while optimizing various objectives of OPF using differential evolution (DE). ∙ Applying ensemble method of CH techniques (ECHT) with SF and SP to handle constraints while optimizing various objectives of OPF using DE. ∙ Statistically compare the results and analyze the performance of SF-DE, SP-DE and ECHT-DE. ∙ Comparison of results of current study with most recent studies on OPF with similar experimental set-up. ∙ Critical analysis of results especially against constraint violation. The organization of rest of the paper is done in following way. Section 2 includes a review of mathematical model including applicable constraints pertaining to OPF in the network. In Section 3, objectives and case studies performed for all the bus test systems are explained with useful numerical values. Description and application of constraint handling (CH) techniques are elaborated in Section 4. Section 5 discusses the simulation results and comparison followed by concluding remarks in Section 6. 82
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
2. Optimal power flow mathematical formulation
2.3.2. Inequality constraints The inequality constraints in OPF reflect the operating limits of the equipment present in the power system and also the limits imposed on lines and load buses to guarantee system security. (a) Generator constraints:
The OPF is a non-linear, non-convex optimization problem and it minimizes certain objectives in the power system subject to many equality and inequality constraints. Mathematically, OPF is represented as:
subject to ∶
𝑉𝐺𝑚𝑖𝑛 ≤ 𝑉𝐺𝑖 ≤ 𝑉𝐺𝑚𝑎𝑥 ∀ 𝑖 ∈ 𝑁𝐺
(6)
𝑃𝐺𝑚𝑖𝑛 ≤ 𝑃𝐺𝑖 ≤ 𝑃𝐺𝑚𝑎𝑥 ∀ 𝑖 ∈ 𝑁𝐺
(7)
𝑚𝑎𝑥 𝑄𝑚𝑖𝑛 ∀ 𝑖 ∈ 𝑁𝐺. 𝐺 ≤ 𝑄𝐺𝑖 ≤ 𝑄𝐺
(8)
𝑖
(1)
Minimize ∶ 𝑓 (𝑥, 𝑢)
𝑖
𝑖
𝑔(𝑥, 𝑢) ≤ 0 ℎ(𝑥, 𝑢) = 0
𝑖
𝑖
𝑖
(b) Transformer constraints:
where, 𝑢 is the vector of control or independent variables, 𝑥 is the vector of state or dependent variables. 𝑓 (𝑥, 𝑢): objective functions of OPF, g (𝑥, 𝑢): set of inequality constraints, ℎ (𝑥, 𝑢): set of equality constraints.
𝑇𝑗𝑚𝑖𝑛 ≤ 𝑇𝑗 ≤ 𝑇𝑗𝑚𝑎𝑥 ∀ 𝑗 ∈ 𝑁𝑇 .
(9)
(c) Shunt compensator constraints: 2.1. Control (independent) variables 𝑚𝑎𝑥 𝑄𝑚𝑖𝑛 ∀ 𝑘 ∈ 𝑁𝐶. 𝐶 ≤ 𝑄 𝐶𝑘 ≤ 𝑄 𝐶 𝑘
The set of variables that can control the power flow in the network is represented in vector form as: 𝑢 = [𝑃𝐺2 … 𝑃𝐺𝑁𝐺 , 𝑉𝐺1 ..𝑉𝐺𝑁𝐺 , 𝑄𝐶1 … 𝑄𝐶𝑁𝐶 , 𝑇1 … 𝑇𝑁𝑇 ]
(2)
(12)
∀ 𝑞 ∈ 𝑛𝑙.
To assess performance of various constraint handling (CH) techniques, several case studies with single and multi-objectives have been performed for standard IEEE 30, 57 and 118-bus test systems. Table 1 includes a summary of major components (viz. generators, transformers, shunt compensators etc.) installed in the test systems and other relevant useful parameters. In 30-bus and 57-bus systems, bus 1 is the swing bus or slack bus, often termed as 𝑉 𝛿 bus. Bus 69 is the swing bus in 118-bus system. In power flow study, the role of swing bus is to balance active and reactive power in the system satisfying power balance Eqs. (4) and (5). For convenience, voltage magnitude (𝑉 ) of swing bus is considered as 1 p.u. while the voltage angle (𝛿) as 0 degree. All other bus voltages and their angles are expressed w.r.t the swing bus, the values of which are obtained as the output of load flow study. Formulation of study cases for the bus test systems is subsequently described.
(3)
3.1. IEEE 30-bus system
2.3. Constraints
A total of ten study cases studies are performed for the system. First six case studies minimize single objective function of OPF. The remaining cases are for multi-objective optimizations which are converted to single objectives with weight factors as in many previous studies and reproduced herein.
As mentioned before, OPF problem has both equality and inequality constraints which are to be satisfied. The constraints are segregated and provided herein. 2.3.1. Equality constraints In OPF, power balance equations are the equality constraints and those are represented as: ( ) 𝑉𝑗 [𝐺𝑖𝑗 cos 𝛿𝑖𝑗 + 𝐵𝑖𝑗 sin(𝛿𝑖𝑗 )] = 0 ∀ 𝑖 ∈ 𝑁𝐵
(4)
( ) 𝑉𝑗 [𝐺𝑖𝑗 sin 𝛿𝑖𝑗 − 𝐵𝑖𝑗 cos(𝛿𝑖𝑗 )] = 0 ∀ 𝑖 ∈ 𝑁𝐵
(5)
3.1.1. Case 1: Minimization of fuel cost This is the most basic objective function of OPF studied almost in all literatures. The association between fuel cost ($/h) and generated power (MW) is approximately given by the quadratic relationship and hence the objective function to be minimized is described as:
𝑗=1 𝑁𝐵 ∑
(11)
3. Objective functions and study cases
where, 𝑃𝐺1 is the generator active power at slack (or swing) bus, 𝑄𝐺𝑖 is the reactive power of generator connected to bus 𝑖, 𝑉𝐿𝑝 is the bus voltage of 𝑝th load bus (PQ bus) and line loading of 𝑞th line is given by 𝑆𝑙𝑞 . 𝑁𝐿 and 𝑛𝑙 are the number of load buses and transmission lines respectively.
𝑄𝐺𝑖 − 𝑄𝐷𝑖 − 𝑉𝑖
𝑝
𝑆𝑙𝑚𝑎𝑥 𝑞
The control variables among the inequality constraints are self-limiting. The optimization algorithm selects a feasible value for each such variable within the defined range. Effective techniques for handling of inequality constraints pertaining to state or dependent variables are discussed in Section 4.
The state of power system is defined by the state variables which can be expressed by vector 𝑥 as:
𝑁𝐵 ∑
≤ 𝑉𝐿𝑝 ≤ 𝑉𝐿𝑚𝑎𝑥 ∀ 𝑝 ∈ 𝑁𝐿
𝑆𝑙𝑞 ≤
2.2. State (dependent) variables
𝑃𝐺𝑖 − 𝑃𝐷𝑖 − 𝑉𝑖
(10)
(d) Security constraints: 𝑉𝐿𝑚𝑖𝑛 𝑝
where, 𝑃𝐺𝑖 is the 𝑖th bus generator active power (except swing generator). Selection of bus 1 as swing bus is representative only and swing bus can be any one of the generator buses. 𝑉𝐺𝑖 is the voltage magnitude at 𝑖th PV bus (generator bus), 𝑇𝑗 is the 𝑗th branch transformer tap, 𝑄𝐶𝑘 is the shunt compensation at 𝑘th bus. 𝑁𝐺, 𝑁𝐶 and 𝑁𝑇 are the number of generators, shunt VAR compensators and transformers respectively. A control variable can assume any value within its range. In reality, transformer taps are not continuous. However, the tap settings expressed here are in p.u. and absolute value of voltage is not accounted for. Hence for study purpose and to compare with past reported results, all control variables including tap settings are considered continuous for most of the study cases. Discrete steps for transformers and shunt capacitors are accounted only in one special study case.
𝑥 = [𝑃𝐺1 , 𝑉𝐿1 … .𝑉𝐿𝑁𝐿 , 𝑄𝐺1 … 𝑄𝐺𝑁𝐺 , 𝑆𝑙1 … 𝑆𝑙𝑛𝑙 ]
𝑘
𝑗=1
𝑓 (𝑥, 𝑢) =
𝑁𝐺 ∑ 𝑖=1
where 𝛿𝑖𝑗 = 𝛿𝑖 − 𝛿𝑗 , is the difference in voltage angles between bus 𝑖 and bus 𝑗, 𝑁𝐵 is the number of buses, 𝑃𝐷 and 𝑄𝐷 are active and reactive load demands, respectively. 𝐺𝑖𝑗 is the transfer conductance and 𝐵𝑖𝑗 is the susceptance between bus 𝑖 and bus 𝑗, respectively.
𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2
𝑖
(13)
where 𝑎𝑖 , 𝑏𝑖 , 𝑐𝑖 are the cost coefficients of the 𝑖th generator producing output power 𝑃𝐺𝑖 . Cost coefficients of IEEE 30-bus system generators are listed in Table A.2 in Appendix. 83
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 1 Summary of systems under study. Items
Buses Branches Generators
IEEE 30-bus system
IEEE 57-bus system
IEEE 118-bus system
Quantity
Details
Quantity
Details
Quantity
Details
30 41 6
Alsac and Stott (1974) Table A.1 (Appendix) Buses: 1 (swing), 2, 5, 8, 11 and 13
57 80 7
Zimmerman et al. (0000) Zimmerman et al. (0000) Buses: 1 (swing), 2, 3, 6, 8, 9 and 12
118 186 54
Zimmerman et al. (0000) Zimmerman et al. (0000) Buses: 1, 4, 6, 8, 10, 12, 15, 18, 19, 24, 25, 26, 27, 31, 32, 34, 36, 40, 42, 46, 49, 54, 55, 56, 59, 61, 62, 65, 66, 69 (swing), 70, 72, 73, 74, 76, 77, 80, 85, 87, 89, 90, 91, 92, 99, 100, 103, 104, 105, 107, 110, 111, 112, 113 and 116 Buses: 5, 34, 37, 44, 45, 46, 48, 74, 79, 82, 83, 105, 107 and 110
Shunt VAR compensation
9
Transformer with tap changer
4
Control variables Connected load Load bus voltage range allowed
24 – 24
Buses: 10, 12, 15, 17, 20, 21, 23, 24 and 29 Branches: 11, 12, 15 and 36
– 283.4 MW, 126.2 MVAr [0.95 – 1.05] p.u.
3
Buses: 18, 25 and 53
17
Branches: 19, 20, 31, 35, 36, 37, 41, 46, 54, 58, 59, 65, 66, 71, 73, 76 and 80 – 1250.8 MW, 336.4 MVAr [0.94 – 1.06] p.u.
33 – 50
14
9
130 – 64
Branches: 8, 32, 36, 51, 93, 95, 102, 107 and 127
– 4242 MW, 1439 MVAr [0.94 – 1.06] p.u.
3.1.2. Case 2: Minimization of cost considering multi-fuels Thermal generating plants have multi-fuel sources such as coal, natural gas and oil for different power output ranges of a unit. In such cases, the fuel cost function or curve is divided into piecewise quadratic cost functions (shown in Fig. 1) depending on the number and nature of fuels used. The cost considering multi-fuels for 𝑖th generator is mathematically expressed as: 𝑓𝑖 (𝑥, 𝑢) = 𝑎𝑖𝑘 + 𝑏𝑖𝑘 𝑃𝐺𝑖 + 𝑐𝑖𝑘 𝑃𝐺2 for 𝑓 𝑢𝑒𝑙𝑘 𝑖
(14)
within power output range 𝑃𝐺𝑚𝑖𝑛 ≤ 𝑃𝐺𝑖 ≤ 𝑃𝐺𝑚𝑎𝑥 ; 𝑘 being the fuel option. 𝑖𝑘 𝑖𝑘 The objective function reflecting the total cost is given by: ) (𝑁𝐺 ∑ (15) 𝑓𝑖 (𝑥, 𝑢) . 𝑓 (𝑥, 𝑢) = 𝑖=1
For the system under study, first two generating units are considered to have multi-fuel options. The corresponding coefficients and power output ranges (MW) are given in Table A.3 in Appendix. Costs for remaining 4-units are same as in case 1. 3.1.3. Case 3: Enhancement of voltage stability of the network Voltage stability problems are receiving growing attention in power systems as network collapses have been experienced in past due to voltage instability. Under normal condition and after being subjected to disturbance, the stability of a power system is characterized by its ability to maintain all bus voltages within acceptable limits. A system enters into a state of voltage instability when a disturbance, increase in load demand or change in system condition causes a progressive and uncontrollable decrease in voltage (Bhattacharya and Chattopadhyay, 2011). Systems with long transmission lines and heavy loading are more prone to voltage instability problem. In power system, improving voltage stability of a system is an important aspect. 𝐿-index of each bus serves as a good indicator of power system stability (Kessel and Glavitsch, 1986). The value of the index varies from 0 to 1, with 0 being the no load case while 1 signifies voltage collapse. If a power system has NL number of load (PQ) buses and NG number of generator (PV) buses, the value of L-index 𝐿𝑗 of bus 𝑗 is defined as: | | 𝑁𝐺 ∑ | 𝑉 | 𝐿𝑗 = ||1 − 𝐹𝑗𝑖 𝑖 || , where 𝑗 = 1, 2, … , 𝑁𝐿 𝑉𝑗 | | 𝑖=1 | | [ ]−1 and 𝐹𝑗𝑖 = − 𝑌𝐿𝐿 [𝑌𝐿𝐺 ]
Fig. 1. Comparison of cost curves of a generator with single fuel and with multi-fuels.
in (17). ][ ] [ ] [ 𝐼𝐿 𝑌 𝑌 𝑉𝐿 = 𝐿𝐿 𝐿𝐺 . 𝐼𝐺 𝑌𝐺𝐿 𝑌𝐺𝐿 𝑉𝐺
(17)
The L-index is calculated for all load buses and maximum value out of those acts as the global indicator for the system stability. Therefore, the objective function of system stability is given by: (18)
𝑓 (𝑥, 𝑢) = 𝐿𝑚𝑎𝑥 = max(𝐿𝑗 ), 𝑤ℎ𝑒𝑟𝑒 𝑗 = 1, 2, …, 𝑁𝐿.
3.1.4. Case 4: Minimization of emission Generation of electrical power from conventional sources of energy emits harmful gases into the environment. The amount of SOx , NOx emission in tonnes per hr (t/h) increases with increase in generated power (in p.u. MW) following the relationship given in Eq. (19). Minimization of emission is set as the objective of OPF.
(16)
where, sub-matrices 𝑌𝐿𝐿 and 𝑌𝐿𝐺 are obtained from system YBUS matrix after separating load (PQ) buses and generator (PV) buses as represented
𝑓 (𝑥, 𝑢) = Emission =
𝑁𝐺 [ ∑ 𝑖=1
84
(𝛼𝑖 + 𝛽𝑖 𝑃𝐺𝑖 + 𝛾𝑖 𝑃𝐺2 ) × 0.01 + 𝜔𝑖 𝑒(𝜇𝑖 𝑃𝐺𝑖 ) 𝑖
]
(19)
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
3.1.8. Case 8: Minimization of fuel cost and voltage deviation Voltage deviation is a measure of voltage quality in the network. The index of deviation is also important from security aspect. The indicator is formulated as cumulative deviation of voltages of all load buses (PQ buses) in the network from nominal value of unity. Mathematically it is expressed as: (𝑁𝐿 ) ∑| | 𝑉𝐷 = (23) |𝑉𝐿𝑝 − 1| . | | 𝑝=1
The combined objective function of fuel cost and voltage deviation is: 𝑓 (𝑥, 𝑢) =
(𝑁𝐺 ∑ 𝑖=1
) 𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2 𝑖
+ 𝜆𝑉 𝐷 × 𝑉 𝐷
(24)
where, weight factor 𝜆𝑉 𝐷 is assigned a value of 100 as in Bouchekara et al. (2016b) and Mohamed et al. (2017). 3.1.9. Case 9: Minimization of fuel cost and enhancement of voltage stability This objective function aims to minimize both fuel cost and enhance voltage stability of the system. The multiple objectives are converted to single objective as: ) (𝑁𝐺 ∑ 𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2 + 𝜆𝐿 × 𝐿𝑚𝑎𝑥 (25) 𝑓 (𝑥, 𝑢) =
Fig. 2. Cost curves of a generator with and without valve point effect.
𝑖=1
𝑖
where, 𝛼𝑖 , 𝛽𝑖 , 𝛾𝑖 , 𝜔𝑖 and 𝜇𝑖 are all emission coefficients provided in Table A.2 in Appendix for 30-bus system.
where, 𝐿𝑚𝑎𝑥 is calculated using Eq. (18) and selected value of weight factor 𝜆𝐿 is 100 (Mohamed et al., 2017).
3.1.5. Case 5: Minimization of real power loss The power loss in transmission system is unavoidable as the lines have inherent resistance. The real power loss (in MW) to be minimized is expressed as:
3.1.10. Case 10: Minimization of fuel cost, emission, voltage deviation and losses Four objectives are combined in this case study. Fuel cost, emission, voltage deviation and real power loss in the network are all minimized simultaneously. The objective function is given by: ) (𝑁𝐺 ∑ 𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2 𝑓 (𝑥, 𝑢) =
𝑓 (𝑥, 𝑢) = 𝑃𝑙𝑜𝑠𝑠 =
𝑛𝑙 ∑
[ ] 𝐺𝑞(𝑖𝑗) 𝑉𝑖2 + 𝑉𝑗2 − 2𝑉𝑖 𝑉𝑗 cos(𝛿𝑖𝑗 )
(20)
𝑞=1
where, 𝛿𝑖𝑗 = 𝛿𝑖 − 𝛿𝑗 , is the difference in voltage angles between bus 𝑖 and bus 𝑗 and 𝐺𝑞(𝑖𝑗) is the transfer conductance of branch 𝑞 connecting buses 𝑖 and 𝑗.
𝑖=1
( ( ))| | 𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2 + ||𝑑𝑖 × sin 𝑒𝑖 × 𝑃𝐺𝑚𝑖𝑛 − 𝑃𝐺𝑖 || 𝑖 𝑖 | | 𝑖=1
𝑁𝐺 ∑
3.2. IEEE 57-bus system General system data of 57-bus system are provided in Table 1. A total of four case studies are performed for the system with half of those being single objective and rest half are multi-objective optimizations. 3.2.1. Case 11: Minimization of fuel cost This basic objective function aims to minimize fuel cost for 57-bus system with optimal power flow. The equation of objective is same as (13). The values of generator cost (Zimmerman et al., 0000) and emission coefficients (Chaib et al., 2016), used in the calculation of cost and emission respectively, are provided in Table A.4 in Appendix.
(21)
where, 𝑑𝑖 and 𝑒𝑖 are the coefficients that represent the valve-point loading effect. The factors used for calculations are provided in Table A.2 in Appendix.
3.2.2. Case 12: Minimization of fuel cost and voltage deviation The goal of the objective function is to minimize both fuel cost and voltage deviation simultaneously. The converted single objective function follows Eq. (24) with weight factor 𝜆𝑉 𝐷 is assigned a value of 100, same as in the case for 30-bus system.
3.1.7. Case 7: Minimization of fuel cost and real power loss The multi-objective case of minimizing fuel cost and real power loss is converted to single objective by multiplying a weight factor to one of the objectives. It is represented as: 𝑓 (𝑥, 𝑢) =
𝑁𝐺 ∑ 𝑖=1
𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2 + 𝜆𝑃 × 𝑃𝑙𝑜𝑠𝑠 𝑖
(26)
The weight factors are selected as in Mohamed et al. (2017) with 𝜆𝐸 = 19, 𝜆𝑉 𝐷 = 21 and 𝜆𝑃 = 22 to balance among the objectives.
3.1.6. Case 6: Minimization of fuel cost considering valve point effect Valve-point effect needs to be considered for more realistic and precise modeling of fuel cost function. The generating units with multi-valve steam turbines exhibit a greater variation in the fuel-cost functions (Bouchekara et al., 2016b). The valve loading effect of multivalve steam turbines is modeled as sinusoidal function, the absolute value of which is added to the basic cost function. The real cost curve function of the steam plant becomes non-continuous as in Fig. 2. The objective of minimization of generating fuel cost with valvepoint effect is given by Biswas et al. (2017): 𝑓 (𝑥, 𝑢) =
𝑖
+ 𝜆𝐸 × 𝐸𝑚𝑖𝑠𝑠𝑖𝑜𝑛 + 𝜆𝑉 𝐷 × 𝑉 𝐷 + 𝜆𝑃 × 𝑃𝑙𝑜𝑠𝑠 .
3.2.3. Case 13: Minimization of fuel cost and enhancement of voltage stability The formulation of the objective function, consisting of both fuel cost and voltage stability, is same as in Eq. (25). The selected weight factor 𝜆𝐿 is also 100.
(22)
where, 𝑃𝑙𝑜𝑠𝑠 is the real power loss in the network calculated by using Eq. (20) and value of factor 𝜆𝑃 is chosen as 40, same as in Mohamed et al. (2017). 85
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
∙ 𝑋𝑖 is feasible but 𝑋𝑗 is infeasible ∙ 𝑋𝑖 and 𝑋𝑗 are both feasible and 𝑋𝑖 has better fitness value than 𝑋𝑗 ∙ 𝑋𝑖 and 𝑋𝑗 are both infeasible, but 𝑋𝑖 has a smaller overall constraint violation (computed by Eq. (29))
3.2.4. Case 14: Minimization of voltage deviation The objective function of minimizing cumulative voltage deviation of load buses from nominal value of 1 p.u. follows Eq. (23). 3.3. IEEE 118-bus system
Therefore, feasible individuals are always considered better than the infeasible individuals in this technique. Two feasible solutions are compared based on their objective function values only while two infeasible solutions are compared based on their overall constraint violations only. Comparison of infeasible solutions based on the overall constraint violation targets to push the infeasible solutions towards feasible region. Comparison of two feasible solutions on the objective value facilitates the convergence to optimal solution.
General system data of 118-bus system are provided in Table 1. A couple of single objective study cases are performed for the system. 3.3.1. Case 15: Minimization of fuel cost This basic objective function of minimizing fuel cost follows Eq. (13). The values of generator cost coefficients are derived from Zimmerman et al. (0000). 3.3.2. Case 16: Minimization of real power loss The objective aims to minimize real power loss in the 118-bus network. Power loss is calculated using Eq. (20). Network data are all taken from Zimmerman et al. (0000).
(b) Self-adaptive penalty (SP) Self-adaptive penalty (Farmani and Wright, 2003; Tessema and Yen, 2006) is an improved version of static penalty. In self-adaptive penalty, two types of penalties are added to each infeasible individual to identify the best infeasible individuals in the current population. The magnitudes of penalties added are controlled by the number of feasible individuals currently present in the combined population. If the number of feasible individuals is less, a higher amount of penalty is added to infeasible individuals with a higher amount of constraint violation. On contrary, if the number of feasible individuals is significantly large, then infeasible individuals with high fitness values will have small penalties added to their fitness values. These two penalties will allow the algorithm to switch between finding more feasible solutions and searching for the optimum solution at any time during the search process. Unlike static penalty, self-adaptive penalty requires no parameter tuning. The final fitness value based on which the population members are ranked is given as 𝐹 (𝑋) = 𝑑 (𝑋) + 𝑝(𝑋), where 𝑑 (𝑋) is the distance value and 𝑝 (𝑋) is the penalty value. The distance value is computed as follows: { 𝜈 (𝑋) , if 𝑟𝑓 = 0 (30) 𝑑 (𝑋) = √ 2 ′′ 2 [𝑓 (𝑋)] + [𝜈 (𝑋)] , otherwise 𝑓 (𝑋) − 𝑓𝑚𝑖𝑛 𝑓 ′′ (𝑋) = (31) 𝑓𝑚𝑎𝑥 − 𝑓𝑚𝑖𝑛 Number of feasible individuals (32) 𝑟𝑓 = Population size
4. Various constraint handling (CH) techniques integrated with DE As evolutionary algorithms (EAs) are stochastic in nature, infeasible solutions should not be straightaway discarded. A proper constraint handling method when used in conjunction with compatible evolutionary algorithm, it can drive the search process towards global feasible optima by making use of the information present in the infeasible solutions. Among the different constraint handling methods, the most commonly employed method is static penalty, where a penalty is added to the fitness of an infeasible solution for violating the constraints. Despite of its simplicity and ease of implementation, the performance of static penalty method heavily depends on the penalty factor which needs to be fine-tuned by trial and error. Therefore, few constraint handling methods that can be employed with EAs have been proposed and are summarized below. A constrained optimization problem with 𝑛 parameters to be optimized is usually written as a nonlinear programming problem of the following form (Qin et al., 2009; Mallipeddi and Suganthan, 2010): ( ) Minimize ∶ 𝑓 (𝑋) , 𝑋 = 𝑥1 , 𝑥2 , … , 𝑥𝑛 and 𝑋 ∈ 𝑆 subject to ∶
(27)
𝑔𝑖 (𝑋) ≤ 0, 𝑖 = 1, … , 𝑝 ℎ𝑗 (𝑋) = 0, 𝑗 = 𝑝 + 1, … , 𝑚
where, 𝜈 (𝑋) is the overall constrain violation as defined in Eq. (29), 𝑓𝑚𝑎𝑥 and 𝑓𝑚𝑖𝑛 are the maximum and minimum values of the objective function 𝑓 (𝑋) in the current combined population. The penalty value is ( ) defined as 𝑝 (𝑋) = 1 − 𝑟𝑓 𝑀 (𝑥) + 𝑟𝑓 𝑁(𝑥) where { 0, if 𝑟𝑓 = 0 𝑀 (𝑋) = and 𝜈 (𝑋) , otherwise { 0, if 𝑋 is a feasible individual 𝑁 (𝑋) = 𝑓 ′′ (𝑋) , if 𝑋is an infeasible individual.
where, 𝑆 is the entire search space, 𝑝 is the number of inequality constraints and (𝑚−𝑝) is the number of equality constraints. The equality constraints are transformed into inequality constraints and thus total constraint is represented as: [ ] { max g𝑖 (𝑋) , 0 , 𝑖 = 1, … , 𝑝 ] [ 𝐺𝑖 (𝑋) = (28) | | max |ℎ𝑖 (𝑋)| − 𝛿, 0 , 𝑖 = 𝑝 + 1, … , 𝑚 where, 𝛿 is a tolerance parameter for the equality constraints. The fitness function 𝑓 (𝑋) is to be minimized so that the optimal solution obtained in the process satisfies the inequality constraint 𝐺𝑖 (𝑋). The overall constraint violation for an infeasible solution is the weighted mean of all constraints and it is defined as: ∑𝑚 𝑤𝑖 (𝐺𝑖 (𝑋)) 𝜈(𝑋) = 𝑖=1∑𝑚 (29) 𝑖=1 𝑤𝑖
(c) Ensemble of constraint handling techniques (ECHT) Generally, each constrained optimization problem is unique in terms of the ratio between feasible search space and the whole search space, multi-modality and the nature of constraint functions. While solving a particular problem, the evolution paths can be different during different runs even if the same algorithm is used because of the stochastic nature of the EAs. In other words, the search process passes through different phases at different points during the evolution. Therefore, depending on several factors such as the ratio between feasible search space and the whole search space, multi-modality of the problem, nature of equality/inequality constraints, the chosen EA and global exploration/local exploitation stages of the search algorithm, different constraint handling methods can be effective during different stages of the search process. Due to the strong interactions among these diverse factors and the stochastic nature of the evolutionary algorithms, it is not straightforward to determine which constraint handling method is the best during a particular stage of the evolution to solve a given
where, 𝑤𝑖 (= 1∕𝐺𝑚𝑎𝑥,𝑖 ) is a weight parameter, 𝐺𝑚𝑎𝑥,𝑖 is the maximum value for violation of constraint 𝐺𝑖 (𝑋) recorded thus far. In the algorithm, 𝑤𝑖 is set as 1∕𝐺𝑚𝑎𝑥,𝑖 , which varies during the evolution to balance the contribution of each constraint in the problem regardless of the different numerical ranges of all constraints. A brief overview of the two constraint handling techniques used for OPF problem is presented herein. (a) Superiority of feasible solutions (SF) In SF (Deb, 2000), solution 𝑋𝑖 is considered superior to solution 𝑋𝑗 when: 86
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
the other constraint handling methods. Therefore, an offspring produced by a particular constraint handling method may be rejected by its own population, but could be accepted by the populations of other constraint handling methods. Hence, in ECHT every function call is utilized effectively making it useful for objective/constraint functions that are computationally expensive. In ECHT, if a particular constraint handling technique is best suited for a given problem and the search method at a particular point in the search process, the offspring population produced by the population of that constraint handling method will dominate the other populations and enter other populations to become parents in the subsequent generations. Therefore, ECHT takes the onus of choosing the best constraint handling technique and tuning the associated parameter values for a particular problem. A summary of steps involved in ECHT-DE is provided in Table 2. Flowchart with basic steps of the algorithm is also presented in Fig. 3. Dimension 𝐷 of OPF problem is 24, 33 and 130 for IEEE 30-bus, 57-bus and 118-bus systems, respectively. Population size 𝑁𝑃 and maximum number of function evaluation, 𝑚𝑎𝑥𝑒𝑣𝑎𝑙 are varied for different cases and those details are mentioned in Section 5. The proposed algorithms are developed using MATLAB software and simulations are carried out on a computer with Intel Core i5 CPU @2.7 GHz and 4GB RAM. Simulation results are discussed in subsequent section.
problem using a given EA. Motivated by these observations, ECHT was developed (Mallipeddi and Suganthan, 2010) to implicitly benefit from the match between constraint handling methods, characteristics of the problem being solved, chosen EA and the exploration–exploitation stages of the search process. (d) Differential evolution (DE) Differential Evolution (DE), introduced by Storn and Price (Storn and Price, 1997), is a stochastic, population based optimization algorithm. The algorithm employs crossover and mutation operators that help the individuals in the population to evolve and improve their fitness. These individuals are evaluated and those that perform better are selected to compose the population in the next generation. A brief description of the steps involved in DE algorithm is presented in this section. Initialization DE optimization process initializes the population of candidate solutions with assignment of random values to each decision vector in the population. These values must be within the feasible bounds (between maximum & minimum) of the decision variables. Initialization of 𝑗th component of 𝑘th decision vector is done as: 𝑋𝑘,0 (𝑗) = 𝑥𝑚𝑖𝑛 (𝑗) + 𝑟𝑎𝑛𝑑𝑘,𝑗 [0, 1] × [𝑥𝑚𝑎𝑥 (𝑗) − 𝑥𝑚𝑖𝑛 (𝑗)]
(33) 5. Results and comparison
where, 𝑟𝑎𝑛𝑑𝑘,𝑗 [0, 1] is a random number lying between 0 and 1; 𝑗 = 1, 2, … , 𝐷 where 𝐷 is the dimension of the decision vector.
Simulation results of application of SF-DE, SP-DE and ECHT-DE algorithms are tabulated, and also detail system by system explanation and comparative analysis are provided in this section. Table 3 summarizes the objectives of various case studies performed for the test systems with all the 3 algorithms. The parameters to be optimized for a specific case are ticked in corresponding cells of the table. A population size (𝑁𝑃 ) of 30 is considered for all the study cases of the systems under study. The population size is selected after several trial runs. Increasing number of population members adversely affects the algorithm in finding feasible solutions. Again, for single objective and multi-objective optimizations of IEEE 30-bus system, 30 000 function evaluations (𝑚𝑎𝑥𝑒𝑣𝑎𝑙 for system 1) are performed in a single complete run of the algorithm. A total of 42 000 function evaluations (𝑚𝑎𝑥𝑒𝑣𝑎𝑙 for system 2) of objective functions are carried out for IEEE 57-bus system in one run. Due to large number of control variables in 118-bus system, 300 000 function evaluations (𝑚𝑎𝑥𝑒𝑣𝑎𝑙 for system 3) are set as the stopping criterion in a study case. Each case of 30-bus and 57-bus systems (i.e. case 1 to case 14) is run 30 times independently to perform statistical analysis and evaluate the performance of the 3-constraint handling techniques.
Mutation In second step, mutation operator of DE creates a donor/mutant vector 𝑉𝑘,𝐺 corresponding to each population member or target vector 𝑋𝑘,𝐺 in the current generation (the subscript ‘G’ denotes parameter at G-th generation). The strategy for mutation used here is ‘DE/rand/1’: 𝑉𝑘,𝐺 = 𝑋𝑅𝑘1,𝐺 + 𝐹 (𝑋𝑅𝑘2,𝐺 − 𝑋𝑅𝑘3,𝐺 ).
(34)
The indices 𝑅𝑘1, 𝑅𝑘2 and 𝑅𝑘3 are mutually exclusive and chosen randomly from the population range. The scaling (mutation) factor 𝐹 is a positive control parameter for scaling the difference vectors. Crossover The donor vector mixes its components with the target vector 𝑋𝑘,𝐺 to form the trial/offspring vector 𝑈𝑘,𝐺 = [𝑈𝑘,𝐺 (1) , 𝑈𝑘,𝐺 (2) , … .., 𝑈𝑘,𝐺 (𝐷)] through crossover. Binomial crossover is adopted here and the scheme operates on each variable whenever a randomly generated number between 0 and 1 is less than or equal to a pre-fixed crossover rate 𝐶𝑅. For an element, the scheme is expressed as: { 𝑉𝑘,𝐺 (𝑗) if 𝑗 = 𝑗𝑟𝑎𝑛𝑑 or 𝑟𝑎𝑛𝑑𝑘,𝑗 [0, 1] ≤ 𝐶𝑅, 𝑈𝑘,𝐺 (𝑗) = (35) 𝑋𝑘,𝐺 (𝑗) otherwise, where 𝑗 = 1, … , 𝐷
5.1. Statistical comparison among ECHT-DE, SF-DE and SP-DE
where, 𝑗𝑟𝑎𝑛𝑑 is a randomly chosen natural number in {1, 2, . . . , D }. The new offspring vector, 𝑈𝑘,𝐺 is evaluated and compared against parent vector 𝑋𝑘,𝐺 in terms of its fitness value and constraint violation. In case of single constraint handling method applied to a problem, the ordering of feasible solutions and operation on infeasible solutions among offspring and parent vectors are dictated by the rules set up by that constraint handling method, i.e. either by SF or SP as described in clause 4(a) and (b). In real-world problem such as OPF where a single run to find optimum settings of control variables with metaheuristics can take considerable time and resource, finding a suitable constraint handling method by trial-and-error may become difficult. The computation time needed to find an appropriate constraint handling method can be saved by using the ECHT method where two constraint handling techniques discussed in Section 4(a) and (b) and DE (Storn and Price, 1997) as the basic search algorithm are used. In ECHT, each constraint handling technique has its own population and parameters that produce offsprings for evaluation. The parent population corresponding to a particular constraint handling method not only competes with its own offspring population but also with offspring populations of
The statistical summary of 30 runs for each study case performed (i.e. case 1 to case 14) using different constraint handling techniques with DE as the basic search algorithm is presented in Table 4. The columns indicate the best, worst, mean and standard deviation values for the objective function in each case. It is evident that no single method is able to provide best fitness or best mean values in all the cases. Further, we carry out Wilcoxon signed rank test to compare performance between any two algorithms. The procedure applied for comparative performance analysis between algorithm ‘A’ and algorithm ‘B’ on a specific study case using signed rank test (Derrac et al., 2011) is described below. ➢ Gather all fitness values over 30 runs of the objective in a case study for both the algorithms ➢ Compute R+, the sum of ranks for runs in which algorithm ‘A’ outperforms algorithm ‘B’ ➢ Compute R−, the sum of ranks for runs in which algorithm ‘B’ outperforms algorithm ‘A’ 87
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 2 Steps in ECHT-DE. ECHT-DE INPUT
∙ Dimension of the problem, 𝐷 ∙ Number of population for each constraint handling method, 𝑁𝑃 ∙ Stopping criteria, maximum number of function evaluation, 𝑚𝑎𝑥𝑒𝑣𝑎𝑙 ∙ Maximum and minimum values of 𝐷-decision variables, in vector form 𝑋𝑚𝑎𝑥 and 𝑋𝑚in 1 D 1 D ∙ 𝑋𝑚𝑎𝑥 = [𝑋𝑚𝑎𝑥 , … , 𝑋𝑚𝑎𝑥 ] and 𝑋𝑚in = [𝑋min , … , 𝑋min ]
INITIALIZATION
∙ Set generation counter, 𝐺 = 0 and function evaluation counter 𝑛𝑓 𝑒𝑣𝑎𝑙 = 0. ∙ Set DE parameters: Mutation factor, 𝐹 = 0.7; Crossover rate, 𝐶𝑅 = 0.9 ∙ POP1: Create population 1 of 𝑁𝑃 individuals uniformly distributed between [𝑋𝑚𝑎𝑥 , 𝑋𝑚𝑖𝑛 ] ∙ POP2: Create population 2 of 𝑁𝑃 individuals uniformly distributed between [𝑋𝑚𝑎𝑥 , 𝑋𝑚𝑖𝑛 ] ∙ Each of the two constraint handling techniques (SF, SP) has its own population POP1 and POP2 ∙ Initialize parameters PAR1 and PAR2 for SF and SP respectively according to Section 4(a) and (b)
EVALUATION
∙ Evaluate objective function, constraint function and constraint violation using Eq. (27) to Eq. (29) for each individual 𝑋1k , ∀𝑘 ∈ {1, … , 𝑁𝑃 } of POP1 ∙ Evaluate objective function, constraint function and constraint violation using Eq. (27) to Eq. (29) for each individual 𝑋2k , ∀𝑘 ∈ {1, … , 𝑁𝑃 } of POP2 ∙ Increase function evaluation counter 𝑛𝑓 𝑒𝑣𝑎𝑙 by 2 × 𝑁𝑃
ALGORITHM LOOP: STEP 1
∙ Update parameter values PAR1 for SF constraint handling methods following Section 4(a) ∙ Update parameter values PAR2 for SP constraint handling methods following Section 4(b)
STEP 2
∙ Perform mutation Eq. (34) and crossover Eq. (35) for each individual of POP1; generate offspring OFS1 ∙ Perform mutation Eq. (34) and crossover Eq. (35) for each individual of POP2; generate offspring OFS2 ∙ Sample mutation and crossover strategy for POP1 at generation 𝐺: Mutant vector, 𝑉 1𝑘,𝐺 = 𝑋1𝑅𝑘1,𝐺 + 𝐹 (𝑋1𝑅𝑘2,𝐺 − 𝑋1𝑅𝑘3,𝐺 ), 𝑘 ∈ {1, … , 𝑁𝑃 } { 𝑉 1𝑘,𝐺 (𝑗) if 𝑗 = 𝑗𝑟𝑎𝑛𝑑 or 𝑟𝑎𝑛𝑑𝑘,𝑗 [0, 1] ≤ 𝐶𝑅, Trial vector element, 𝑈 1𝑘,𝐺 (𝑗) = 𝑋1𝑘,𝐺 (𝑗) otherwise, where 𝑗 = 1, … , 𝐷 where 𝑅𝑘1, 𝑅𝑘2 and 𝑅𝑘3 are mutually exclusive integers randomly chosen from population range {1, … , 𝑁𝑃 }; 𝑗𝑟𝑎𝑛𝑑 is a randomly chosen natural number within {1, … , 𝐷}
STEP 3
∙ Compute objective function, constraint function and constraint violation using Eq. (27) to Eq. (29) for each individual 𝑋1′k , ∀𝑘 ∈ {1, … , 𝑁𝑃 } of OFS1 ∙ Compute objective function, constraint function and constraint violation using Eq. (27) to Eq. (29) for each individual 𝑋1′k , ∀𝑘 ∈ {1, … , 𝑁𝑃 } of OFS2 ∙ Increase function evaluation counter 𝑛𝑓 𝑒𝑣𝑎𝑙 by 2 × 𝑁𝑃
STEP 4
∙ Combine parent population with all offsprings to create groups. ∙ Group1: (POP1+OFS1+OFS2) and Group2: (POP2+OFS1+OFS2)
STEP 5
∙ In selection step, parent populations POP1 and POP2 for the next generation are selected from Group1 and Group2 respectively. In a Group (say Group1), since OFS1 is produced by POP1 by mutation and crossover, DE’s selection based on competition between parent and its offspring is employed when POP1 competes with OFS1. But when POP1 competes with OFS2, produced by other populations, each member in POP1 competes with a randomly selected offspring from OFS2.
STEP 6
∙ Is stopping criterion 𝑚𝑎𝑥𝑒𝑣𝑎𝑙 reached? ∙ If yes, STOP ∙ If not, increase generation counter by 1, i.e. 𝐺 = 𝐺 + 1. Go to algorithm loop STEP 1.
Table 3 Summary of case studies. Case no.
IEEE 30-bus system Basic fuel cost
Case 1 Case 2 Case 3 Case 4 Case 5 Case 6 Case 7 Case 8 Case 9 Case 10 Case 11 Case 12 Case 13 Case 14 Case 15 Case 16
Multi-fuel cost
IEEE 57-bus system Voltage stability
Emission
Power loss
Voltage deviation
Valve-point effect
Basic fuel cost
Voltage stability
✔ ✔ ✔
✔
IEEE 118-bus system Voltage deviation
Basic fuel cost
Power loss
✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔
✔ ✔ ✔ ✔ ✔
✔
✔ ✔ ✔ ✔ ✔
➢ Calculate 𝑝-value which determines the significance of results in a statistical hypothesis test. The smaller the 𝑝-value, the stronger is the proof against null hypothesis 𝐻0.
8. In the remaining cases, SF-DE is found to be superior. ECHT-DE outperforms SP-DE in cases 1, 3, 6, 7 and 9 for 30-bus system. However, for all cases of 57-bus system, SP-DE outdoes ECHT-DE. An interesting point to note is in the performance comparison between SF-DE and SP-DE. While in cases 1, 3, 6, 7 and 9, SP-DE is inferior to SF-DE, in all remaining cases it is equal or better. This very fact reiterates the need and justification for development of ensemble method. The factors influencing the performance of a CH technique are the nature of the objective function, constraints, the feasible and infeasible regions for the
The results obtained using Wilcoxon signed rank test are presented in Table 5. The column 𝐻0 defines whether null hypothesis is valid or not. If null hypothesis is valid (i.e. 𝐻0 = ‘Yes’ with a significance level, 𝛼 = 0.05), the performance of two methods is statistically same for the study case. ECHT-DE and SF-DE perform equally in cases 4 and 88
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Fig. 3. Flowchart for implementation of ECHT-DE.
Table 4 Statistical summary of case studies using 3-CH techniques. Case no.
ECHT-DE
SF-DE
SP-DE
Best
Worst
Mean
Std dev
Best
Worst
Mean
Std dev
Best
Worst
Mean
Case 1
800.4148
800.4258
800.4206
0.0026
800.4131
800.4192
800.4151
0.0015
800.4293
800.4684
800.4413
0.0100
Case 2
646.4532
646.5423
646.4958
0.0232
646.4390
646.5117
646.4776
0.0178
646.4314
646.496
646.4565
0.0169
Case 3
0.13632
0.13721
0.13693
0.0002
0.13671
0.13704
0.13685
0.0001
0.13743
0.13857
0.13777
0.0002
Case 4
0.20482
0.20482
0.20482
0
0.20482
0.20482
0.20482
0
0.20482
0.20482
0.20482
0
Case 5
3.0850
3.0871
3.0858
0.0005
3.0845
3.0857
3.0849
0.0003
3.0844
3.08536
3.0848
0.0003
Case 6
832.1356
832.2239
832.1811
0.0222
832.0882
832.1291
832.1056
0.0105
832.4813
832.876
832.655
0.0963
Case 7
1040.151
1040.233
1040.181
0.0213
1040.125
1040.162
1040.14
0.0096
1040.134
1040.337
1040.239
0.0444
Case 8
813.1742
813.4095
813.2470
0.049
813.1956
813.3376
813.2585
0.0361
813.1959
813.2643
813.2306
0.0181
Case 9
814.1708
814.2001
814.1843
0.0075
814.1649
814.1956
814.1767
0.0063
814.1841
814.2273
814.2017
0.0121
Case 10
964.1331
964.1564
964.1437
0.0061
964.1254
964.1418
964.1307
0.0038
964.1234
964.1399
964.1276
0.0034
Case 11
41670.56
41680.09
41675.33
2.2866
41667.85
41684.29
41671.16
3.5896
41667.82
41678.96
41670.44
2.2474
Case 12
41776.48
41783.22
41778.87
1.7681
41775.09
41778.39
41776.16
0.6806
41774.75
41776.19
41775.3
0.3724
Case 13
41699.25
41727.03
41706.5
6.1009
41695.55
41712.35
41699.51
3.7249
41696.54
41700.97
41698.17
1.0153
Case 14
0.60416
0.64184
0.62012
0.0094
0.59584
0.61677
0.60602
0.0059
0.59267
0.60284
0.59800
0.0028
solutions, and the way the CH technique works. In SP, higher amount of penalty is added to infeasible individuals if a few feasible individuals are present in the population (Tessema and Yen, 2006). On the other hand, basic SF technique has the problem of maintaining population diversity and may stagnate especially when ratio between feasible region and search space is small and initial population has a few feasible solutions (Coello, 2002). The limited number of feasible solutions during initial phases of search process can deter SF from attaining
Std dev
optimum values in the trial runs. In summary, it is a tedious task to select appropriate CH method for a practical problem evaluated under different scenarios. Gauging extent of non-linearity in the objective function or assessing nature of constraints or delimiting feasible region is not straightforward. Ensemble method has the advantage of accessing and operating on the population members of all CH techniques it is composed of. Due to probabilistic nature of DE, ECHT-DE is able to attain close to, if not the best fitness value in most trial runs. 89
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 5 Wilcoxon signed rank test results. Case no. Case 1 Case 2 Case 3 Case 4 Case 5 Case 6 Case 7 Case 8 Case 9 Case 10 Case 11 Case 12 Case 13 Case 14
ECHT-DE vs SF-DE
ECHT-DE vs SP-DE
SF-DE vs SP-DE
R+
R−
𝑝-value
𝐻0
R+
R−
𝑝-value
𝐻0
R+
R−
𝑝-value
H0
1 94 103.5 0 8.5 0 0 303 31 5 53 0 11 13.5
464 371 361.5 0 456.5 465 465 162 434 460 412 465 454 451.5
1.92× 10−6 4.39 × 10−3 0.00793 1 4.07 × 10−6 1.73 × 10−6 1.73 × 10−6 0.14704 3.41 × 10−5 2.88 × 10−6 2.22 × 10−4 1.73 × 10−6 5.22 × 10−6 6.65 × 10−6
No No No Yes No No No Yes No No No No No No
465 16 465 0 1 465 446 176 425 0 28 0 2 0
0 449 0 0 464 0 19 289 10 465 437 465 463 465
1.73 × 10−6 8.47 × 10−6 1.73 × 10−6 1 1.92 × 10−6 1.73 × 10−6 1.13 × 10−5 0.24519 7.23 × 10−6 1.73 × 10−6 2.60 × 10−5 1.73 × 10−6 2.13 × 10−6 1.73 × 10−6
No No No Yes No No No Yes No No No No No No
465 67 465 0 162.5 465 464 83 465 72 191 14 168 12
0 398 0 0 302.5 0 1 382 0 393 274 451 297 453
1.73 × 10−6 6.64 × 10−4 1.73 × 10−6 1 0.14986 1.73 × 10−6 1.92 × 10−6 2.10 × 10−3 1.73 × 10−6 9.63 × 10−4 0.39333 6.98 × 10−6 0.18462 5.75 × 10−6
No No No Yes Yes No No No No No Yes No Yes No
5.2. System 1: IEEE 30-bus test system In this sub-section initially in Table 6, we list down settings of all control and state (dependent) variables alongwith their allowable ranges for the best fitness value obtained for the objective function in a case study pertaining to 30-bus test system using one of the 3CH techniques (ECHT, SF and SP). Active power of swing generator (𝑃G1 for 30 and 57-bus systems) and reactive power of all generators are state or dependent variables, treated as constraints in the optimization. Values of these variables are listed here to show that the CH techniques duly comply with the limits of these constraining variables. Limits of allowable reactive power for generators are derived from Zimmerman et al. (0000). Last row of the table indicates the approximate time needed by the processor (CPU) for one run of a study case. Table 7 presents comparison of results obtained using the 3-CH techniques in present study with past studies for single objective functions. Table 8 provides similar comparison for multi-objective cases. For case 1 of optimizing basic fuel cost, ECHT-DE and SF-DE algorithms lead to fuel costs of 800.4148 $/h and 800.4131 $/h, respectively satisfying all system constraints, more precisely conforming to the important inequality constraints on generator reactive power, line capacity and load bus voltage. Among all inequality constraints, constraint on load bus voltage is found to be critical as the operating voltages of load buses are often found to be near the limits. Some of the recent literatures recorded better results than what have been achieved by the 3-methods in present study. A careful observation of those results reveals violation of voltage limits [0.95 p.u.−1.05 p.u.], thus rendering the solutions infeasible. Cursory glance at voltage profiles of load buses leads to the presumption that the violation has all likely occurred exceeding the upper limit. This situation of overvoltage is undesirable as it may stress the system and often lead to failure of the connected equipment. For the 30-bus system, a total of 24 load buses (PQ bus) are present. Theoretically, if all these buses operate at the limits, the maximum possible cumulative absolute value of voltage deviation (VD in Eq. (23)) would be 1.2 p.u. (i.e. 24×0.05 p.u.). The reported VD values in some literatures are higher than 1.2 p.u. in many cases. As it turns out, the contentious results are output of static penalty method which is sensitive to selection of penalty coefficient. In comparison Tables 7 and 8, such cases are highlighted with footnotes. In this specific problem of OPF, limits on all state or dependent variables are duly satisfied in our study. In static penalty function approach, the constraints on these variables are violated in many cases, possibly without cognizance of the programmer. Moreover, the superiority of a proper CH technique lies in its ability to find the best result by operating near to the limits but not violating those. Figs. 4 and 5 indicate the load bus voltages for various study cases performed for IEEE 30-bus system. The voltage profile for a study case corresponds to the control variables obtained using the best performing method on that case as listed in Table 6. Clearly voltage limits are all conformed for the best solution.
Fig. 4. IEEE 30-bus system – voltage profiles of load buses for best solutions of case 1 to case 5.
Fig. 5. IEEE 30-bus system – voltage profiles of load buses for best solutions of case 6 to case 10.
In case 2 of minimizing cost with multiple fuels, the present study with SP-DE finds a cost of 646.4314 $/h and turns out to be among the best values when compared with results achieved by other algorithms. In case 3 of voltage stability indicator, i.e. minimization of maximum ‘L-index’ (𝐿𝑚𝑎𝑥 ) of the system load buses, ECHT-DE gives best output of 0.13632, marginally better than the other equivalent algorithms. Objective of emission in tonnes per hour (t/h) and other useful calculated parameters obtained by ECHT-DE, SF-DE and SP-DE in case 4 are almost same with those reported by algorithms MSA (Mohamed 90
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 6 Simulation results of best solutions for system 1: IEEE 30-bus system. Parameters Algorithm 𝑃G2 (MW) 𝑃G5 (MW) 𝑃G8 (MW) 𝑃G11 (MW) 𝑃G13 (MW) 𝑉1 (p.u.) 𝑉2 (p.u.) 𝑉5 (p.u.) 𝑉8 (p.u.) 𝑉11 (p.u.) 𝑉13 (p.u.) 𝑄c10 (MVAr) 𝑄c12 (MVAr) 𝑄c15 (MVAr) 𝑄c17 (MVAr) 𝑄c20 (MVAr) 𝑄c21 (MVAr) 𝑄c23 (MVAr) 𝑄c24 (MVAr) 𝑄c29 (MVAr) 𝑇11 (p.u.) 𝑇12 (p.u.) 𝑇15 (p.u.) 𝑇36 (p.u.) Fuel cost ($/h) Emission (t/h) 𝑃loss (MW) VD (p.u.) L-index (max) 𝑃G1 (MW) 𝑄G1 (MVAr) 𝑄G2 (MVAr) 𝑄G5 (MVAr) 𝑄G8 (MVAr) 𝑄G11 (MVAr) 𝑄G13 (MVAr) CPU time (s)
Min
Max
20 15 10 10 12 0.95 0.95 0.95 0.95 0.95 0.95 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.90 0.90 0.90 0.90
80 50 35 30 40 1.10 1.10 1.10 1.10 1.10 1.10 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 1.10 1.10 1.10 1.10
50 −20 −20 −15 −15 −10 −15
200 150 60 62.5 48.7 40 44.7
Case 1
Case 2
Case 3
Case 4
Case 5
Case 6
Case 7
Case 8
Case 9
Case 10
SF-DE 48.7282 21.3535 21.1986 11.9028 12.0021 1.0832 1.064 1.0333 1.038 1.0808 1.0418 3.6349 2.6015 3.8678 4.9883 4.184 4.9792 3.1586 4.9816 2.3958 1.051 0.9254 0.9649 0.9727 800.4131 0.36652 9.0104 0.92097 0.13786 177.2252 2.9154 19.2342 26.0422 28.0047 21.0932 −6.0223 133.1
SP-DE 54.9987 23.9573 34.9706 18.6602 17.5408 1.0753 1.0612 1.0314 1.0406 1.077 1.051 4.3281 0.2885 3.812 4.1986 4.0378 4.7194 4.0527 4.7473 2.3624 1.0231 0.9575 0.9765 0.9767 646.4314 0.28351 6.7276 0.91253 0.13832 140 −2.4871 17.6727 23.9799 27.182 16.2243 0.9936 122.2
ECHT-DE 79.7793 49.718 34.1924 28.8158 13.733 1.0542 1.051 1.0682 1.0568 1.0997 1.0741 2.9321 0.3526 1.5598 1.2874 0.0082 0.0327 0.0021 0.2923 0.0232 1.0445 0.9001 1.0035 0.9646 917.5916 0.2252 4.5224 0.9110 0.13632 81.684 −19.9831 −19.9361 54.341 48.4862 27.0418 18.6036 136.5
ECHT-DE 67.5933 50.0 35.0 30.0 40.0 1.0622 1.0565 1.0373 1.0439 1.0691 1.0507 1.5203 3.8625 4.2725 4.9037 3.8882 4.9426 3.1486 4.9819 2.5522 1.0677 0.901 0.9905 0.9771 944.3782 0.20482 3.2179 0.89240 0.13844 64.0253 −5.6506 8.0392 21.9063 28.345 19.4601 1.5932 138.2
SP-DE 80.0 50.0 35.0 30.0 40.0 1.0613 1.0572 1.0378 1.0442 1.0685 1.055 3.3113 0.8576 4.1219 4.9901 4.0549 5.0 2.8962 4.9919 2.2836 1.0656 0.9071 0.9908 0.9763 967.5962 0.20726 3.0844 0.90359 0.13832 51.4968 −5.5594 7.2517 21.7516 27.7104 18.6164 4.7497 136.4
SF-DE 44.8982 18.5075 10.023 10.0013 12.0006 1.0846 1.062 1.0285 1.0348 1.084 1.0516 2.9496 3.7232 3.9523 4.8533 3.5919 4.9914 3.6448 5.0 2.2825 1.0068 0.9911 0.9796 0.9773 832.0882 0.43730 10.6387 0.84935 0.13934 198.6081 4.9696 16.9485 24.3502 28.5989 18.1039 1.3342 137.6
SF-DE 55.4313 38.1719 35.0 30.0 26.7077 1.0687 1.0581 1.0348 1.0427 1.0753 1.0507 2.8395 0.8649 4.3812 4.9981 3.9664 4.9924 2.7251 4.9889 2.2225 1.0634 0.9076 0.9811 0.9737 859.1458 0.2289 4.5245 0.92731 0.13785 102.6169 −3.3869 10.9563 23.0135 27.5277 20.7233 1.0117 137.9
ECHT-DE 48.6293 21.6905 22.6953 11.8749 12.0233 1.0401 1.024 1.0144 1.0059 1.0424 0.9881 4.9484 0.5117 4.9733 0.0085 4.9924 4.9014 4.9636 4.9872 2.5552 1.0624 0.9002 0.9393 0.97 803.7198 0.36384 9.8414 0.09454 0.14888 176.3282 −4.3297 14.3453 46.9979 40.3091 21.4086 −14.6404 123.3
SF-DE 48.786 21.212 21.5819 11.8277 12.0002 1.0831 1.0641 1.0329 1.0379 1.0928 1.0455 0.7465 0.301 4.517 4.7206 3.7417 4.9716 2.7587 4.9129 2.0903 1.0265 0.9426 0.9667 0.9697 800.4203 0.36592 8.9985 0.93846 0.13745 176.9908 2.5344 20.1818 25.7862 28.8048 22.6468 −3.2481 130.4
SP-DE 52.528 31.4679 35.0 26.6793 21.0893 1.0723 1.0589 1.0317 1.0403 1.0305 1.0124 3.6781 3.9279 4.0292 4.9399 4.9874 5.0 4.1832 5.0 2.5508 1.0882 0.9529 1.021 1.0043 830.2123 0.25294 5.5857 0.29615 0.14756 122.2219 −1.0989 13.1559 22.963 26.8471 15.8058 0.7012 128.8
et al., 2017) and ARCBBO (Kumar and Premalatha, 2015). Case 5 by SF-DE and SP-DE result in real power losses of 3.0845 MW and 3.0844 MW, respectively, comparable with the results summarized in the table. Valve-point effect is considered in case 6 to arrive at an increase in cost than in case 1 with final value being 832.0882 $/h, obtained by SFDE. In summary, although variation in performance is observed among 3-CH techniques, output results of one or more methods employed in our study are among the best reported results on the problem of OPF. Moreover, the objective of this paper is not to prove superiority on mere numerical results, but to demonstrate strict compliance with the system constraints by the CH techniques. Case 7 to case 10 are for multi-objective cases of OPF for 30bus system. In these case studies, fitness of the combined objective function is the major deciding factor in ranking the results out of various optimization algorithms. For a meaningful comparison, fitness value of other methods is calculated and produced here using the same weight factor of the respective objective function. In multi-objective cases, a change in weight factor e.g. higher weight factor on fuel cost in case 7 will bring down the cost but will increase power loss. Quick review of Table 8 reveals that either one or more algorithms among ECHT-DE, SF-DE and SP-DE achieve best fitness values in all the cases. Though best fitness is reported by SF-DE in case 7, an intermediate value fuel cost, one of the constituting objectives, is achieved. The other objective of real power loss is the least when assessed with other algorithms. In case 8 of minimizing cost and voltage deviation (𝑉 𝐷), 𝑉 𝐷 achieved by ECHT-DE is the lowest among all other comparable methods. Optimum values of both fitness and cost are the lowest in run of case 9 by SF-DE algorithm. Four objectives of cost, emission, real power loss and voltage deviation are simultaneously minimized in case 10. Alongwith fitness value, SP-DE arrives at lowest values of cost and loss when compared with MSA (Mohamed et al., 2017) and MFO (Mohamed et al., 2017).
Fig. 6. Comparative convergence of case 1 for IEEE 30-bus system.
Graphical comparisons of convergences of 3-CH techniques for case 1, case 2 and case 6 of fuel cost related objective functions are illustrated in Figs. 6, 7 and 11 respectively. The convergence speeds are not remarkably different among the algorithms. However, rapid and sudden convergence is observed for both ECHT-DE and SP-DE during initial phase of the search process. SF-DE converges to the optimal solution more steadily. For cases 3, 4 and 5, the convergences are indicated in Figs. 8–10 respectively. The convergence is a little more uneven in case 3 due to the nature of the objective function. In general, the algorithms converge to optimal solutions somewhere around 15 000 function evaluations. Convergences of two-objective cases are provided in Fig. 12–14. For clarity, convergence of only one algorithm attaining best fitness 91
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 7 Comparison of 3-CH techniques with previous studies for single objective cases of IEEE 30-bus system. Case #
Algorithm
Fuel cost ($/h)
Emission (t/h)
𝑃loss (MW)
VD (p.u.)
L-index (max)
Case 1
ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) AGSO (Daryani et al., 2016) ARCBBO (Kumar and Premalatha, 2015) SKH (Pulluri et al., 2016) BSA (Chaib et al., 2016) ICBO (Bouchekara et al., 2016b) APFPA (Mahdad and Srairi, 2016) FHSA (Pandiarajan and Babulal, 2016) GEM (Bouchekara et al., 2016a) DE (Shaheen et al., 2016) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) LTLBO (Ghasemi et al., 2015) GABC (Roy and Jadhav, 2015) BSA (Chaib et al., 2016) ICBO (Bouchekara et al., 2016b) ECHT-DE SF-DE SP-DE SKH (Pulluri et al., 2016) GEM (Bouchekara et al., 2016a) DE (Shaheen et al., 2016) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) DSA (Abaci and Yamacli, 2016) AGSO (Daryani et al., 2016) ARCBBO (Kumar and Premalatha, 2015) GEM (Bouchekara et al., 2016a) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) DSA (Abaci and Yamacli, 2016) ARCBBO (Kumar and Premalatha, 2015) APFPA (Mahdad and Srairi, 2016) GEM (Bouchekara et al., 2016a) ECHT-DE SF-DE SP-DE PSO (Bouchekara et al., 2016b) ICBO (Bouchekara et al., 2016b) BSA (Chaib et al., 2016) APFPA (Mahdad and Srairi, 2016)
800.4148 800.4131 800.4293 800.5099 801.75 800.5159 800.5141 799.0760a 799.0353a 798.9144a 799.914a 799.0463a 799.0827a 646.4532 646.4390 646.4314 646.8364 647.4315 647.03 646.1504a 645.1668a 917.5916 875.8929 842.015 814.0100 816.9095a 915.2172a 944.3782 944.3242 944.3323 944.5003 944.4086 953.629 945.1597 943.6358a 967.6001 967.6198 967.5962 967.6636 967.6493 967.6605 965.6590a 966.7473a 832.1356 832.0882 832.4813 832.6871 830.4531a 830.7779a 830.4065a
0.36593 0.36652 0.36806 0.36645 0.3703 0.3663 0.3662 0.3671 – – – 0.3665 – 0.28351 0.2835 0.28351 0.28352 0.2835 – 0.2833 – 0.2252 0.22801 0.24635 0.3740 0.2802 – 0.20482 0.20482 0.20482 0.20482 0.20583 0.2059 0.2048 0.2048a 0.20726 0.20726 0.20726 0.20727 0.20826 0.2073 – 0.2072 0.43765 0.43730 0.43651 – – 0.4377 –
8.9999 9.0104 9.0583 9.0345 – 9.0255 9.0282 8.6543 8.6132 8.5800 – 8.6257 8.63 6.722 6.7258 6.7276 6.8001 6.9347 6.8160 6.6233 6.3828 4.5224 4.6412 5.92 9.9056 6.2313 3.626 3.2179 3.2179 3.2178 3.2358 3.2437 – 3.2624 3.0160 3.0850 3.0845 3.0844 3.1005 3.0945 3.1009 2.8463a 2.8863a 10.6772 10.6387 10.6762 – 10.2370 10.2908 10.2178
0.91231 0.92097 0.91015 0.90357 – 0.8867 – 1.9129a 1.9652a 1.9451a 1.5265a 1.9312a 1.8505a 0.92671 0.93062 0.91253 0.84479 0.8896 0.8010 1.0273a 1.8232a 0.9110 0.90387 0.80674 – 1.8320a 2.1064a 0.89240 0.89617 0.90215 0.87393 – – 0.8647 1.9504a 0.90937 0.90648 0.90359 0.88868 – 0.8913 2.0720a 1.9755a 0.80326 0.84935 0.75042 – 1.7450a 1.2050a 1.8909a
0.13777 0.13786 0.13775 0.13833 – 0.1385 0.1382 0.1273 0.1261 – – 0.1264 0.1277 0.13812 0.13776 0.13832 0.13867 – – 0.1378 0.1282 0.13632 0.13671 0.13743 0.1366 0.1257 0.1243 0.13844 0.13844 0.13830 0.13888 0.12734 – 0.1387 0.1269 0.13836 0.13818 0.13832 0.13858 0.12604 0.1386 – 0.1265 0.14039 0.13934 0.14157 – 0.1289 0.1363 –
Case 2
Case 3
Case 4
Case 5
Case 6
a
Infeasible solution, constraint on load bus voltage is violated.
Fig. 8. Comparative convergence of case 3 for IEEE 30-bus system.
Fig. 7. Comparative convergence of case 2 for IEEE 30-bus system.
92
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 8 Comparison of 3-CH techniques with previous studies for multi-objective cases of IEEE 30-bus system. Case #
Algorithm
Fitness
Fuel cost ($/h)
Emission (t/h)
𝑃loss (MW)
VD (p.u.)
L-index (max)
Case 7
ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) MFO (Mohamed et al., 2017) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) MFO (Mohamed et al., 2017) BSA (Chaib et al., 2016) ICBO (Bouchekara et al., 2016b) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) FPA (Mohamed et al., 2017) BSA (Chaib et al., 2016) ICBO (Bouchekara et al., 2016b) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) MFO (Mohamed et al., 2017)
1040.151 1040.125 1040.134 1040.808 1041.671 813.1742 813.1956 813.1959 814.1545 814.3541 814.8994 813.5378 814.1708 814.1649 814.1841 814.9378 814.9067 812.9240a 811.8477a 964.1331 964.1254 964.1234 965.2905 965.8077
858.867 859.1458 858.9319 859.1915 858.5812 803.7198 803.4241 803.4196 803.3125 803.7911 803.4294 803.3978 800.4321 800.4203 800.4365 801.2248 801.1487 800.3340a 799.3277a 830.1156 830.1366 830.2123 830.639 830.9135
0.22902 0.2289 0.22895 0.22899 0.22947 0.36384 0.36424 0.36324 0.36344 0.36355 0.3546 – 0.36585 0.36592 0.36517 0.36106 0.3718 0.3514 – 0.25293 0.25313 0.25294 0.25258 0.25231
4.5321 4.5245 4.5301 4.5404 4.5772 9.8414 9.7807 9.7573 9.7206 9.8685 9.3751 9.7453 9.0043 8.9985 8.9838 8.9761 9.3174 8.4904 8.6465 5.5894 5.5887 5.5857 5.6219 5.5971
0.93028 0.92731 0.92626 0.92852 0.89944 0.09454 0.09772 0.09776 0.10842 0.10563 0.1147 0.1014 0.91244 0.93846 0.93706 0.92655 0.87563 1.9855a 1.9961a 0.29738 0.29653 0.29615 0.29385 0.33164
0.13796 0.13785 0.13781 0.13814 0.13806 0.14888 0.14893 0.14893 0.14783 0.14906 0.14840 0.14900 0.13739 0.13745 0.13748 0.13713 0.13758 0.12590 0.12520 0.14748 0.14756 0.14756 0.14802 0.14556
Case 8
Case 9
Case 10
a
Infeasible solution, constraint on load bus voltage is violated.
Fig. 11. Comparative convergence of case 6 for IEEE 30-bus system.
Fig. 9. Comparative convergence of case 4 for IEEE 30-bus system.
Fig. 12. Convergence of case 7 (SF-DE) for IEEE 30-bus system. Fig. 10. Comparative convergence of case 5 for IEEE 30-bus system.
5.3. System 2: IEEE 57-bus test system value is included in the diagram. As an obvious fact, fitness function Settings of control and state variables alongwith their ranges for the best fitness values of study cases for 57-bus test system are listed in Table 9. In Table 10, the results of all 3-CH techniques (ECHT, SF and SP) are compared with most recent studies available for same objective functions. As it turns out, SP-DE leads to best fitness values in all the
gradually reduces in all these cases, though the individual objectives show irregularities during the process of convergence. The fact is due to non-linear relationship of the individual objective functions with control parameters of the network. 93
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 9 Simulation results of best solutions for system 2: IEEE 57-bus system. Parameters Algorithm 𝑃G2 (MW) 𝑃G3 (MW) 𝑃G6 (MW) 𝑃G8 (MW) 𝑃G9 (MW) 𝑃G12 (MW) 𝑉1 (p.u.) 𝑉2 (p.u.) 𝑉3 (p.u.) 𝑉6 (p.u.) 𝑉8 (p.u.) 𝑉9 (p.u.) 𝑉12 (p.u.) 𝑄c18 (MVAr) 𝑄c25 (MVAr) 𝑄c53 (MVAr) 𝑇19 (p.u.) 𝑇20 (p.u.) 𝑇31 (p.u.) 𝑇35 (p.u.) 𝑇36 (p.u.) 𝑇37 (p.u.) 𝑇41 (p.u.)
Min 30 40 30 100 30 100 0.95 0.95 0.95 0.95 0.95 0.95 0.95 0 0 0 0.90 0.90 0.90 0.90 0.90 0.90 0.90
Max
Case 11
Case 12
Case 13
Case 14
Parameters
Min
Max
Case 11
Case 12
Case 13
Case 14
100 140 100 550 100 410 1.10 1.10 1.10 1.10 1.10 1.10 1.10 20 20 20 1.10 1.10 1.10 1.10 1.10 1.10 1.10
SP-DE 89.759 44.564 71.5502 459.8236 98.2124 359.4297 1.0645 1.062 1.0548 1.0598 1.0754 1.0488 1.0481 7.556 13.2503 11.4127 1.046 0.9656 1.0192 1.0023 1.0356 1.0294 0.9955
SP-DE 89.1352 45.135 71.6172 460.5427 96.338 360.8131 1.0328 1.0311 1.0262 1.0414 1.0624 1.0277 1.0179 4.9653 14.9616 14.9965 0.9982 0.9796 0.9918 1.0358 0.9643 1.0239 1.0151
SF-DE 90.2581 44.907 71.1695 461.2539 94.1207 360.8523 1.062 1.0598 1.0528 1.0591 1.0752 1.0477 1.0473 7.5494 13.3646 12.0559 0.9971 0.9803 1.0104 1.0088 1.0149 1.0315 0.9964
SP-DE 68.1397 92.9072 34.0623 294.8929 95.3897 373.4152 1.0086 1.0065 1.0132 1.0021 1.0176 1.0082 1.0392 4.0015 19.7885 19.9272 0.923 1.0505 0.9712 1.0577 1.0624 1.0043 0.9926
𝑇46 (p.u.) 𝑇54 (p.u.) 𝑇58 (p.u.) 𝑇59 (p.u.) 𝑇65 (p.u.) 𝑇66 (p.u.) 𝑇71 (p.u.) 𝑇73 (p.u.) 𝑇76 (p.u.) 𝑇80 (p.u.) Fuel cost ($/h) Emission (t/hr) 𝑃loss (MW) VD (p.u.) L-index (max) 𝑃G1 (MW) 𝑄G1 (MVAr) 𝑄G2 (MVAr) 𝑄G3 (MVAr) 𝑄G6 (MVAr) 𝑄G8 (MVAr) 𝑄G9 (MVAr) 𝑄12 (MVAr) CPU time (s)
0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.90
1.10 1.10 1.10 1.10 1.10 1.10 1.10 1.10 1.10 1.10
0 −140 −17 −10 −8 −140 −3 −150
576 200 50 60 25 200 9 155
0.9543 0.9134 0.9783 0.9696 0.9754 0.9421 0.9743 0.9952 0.9595 0.9987 41667.82 1.350 14.9090 1.54367 0.28123 142.37 47.7018 48.6597 35.7354 −6.5273 54.45 8.7382 53.0263 219.9
0.9403 0.9015 0.966 0.9667 0.982 0.9349 0.9666 1.0007 0.9411 1.0075 41697.50 1.3550 15.5897 0.77253 0.29228 142.8084 42.0835 49.8259 33.8151 −4.9521 75.3832 8.8451 44.529 203.6
0.9638 0.9127 0.9772 0.9651 0.975 0.9397 0.9734 0.9911 0.9637 0.9921 41667.53 1.3576 14.8963 1.61174 0.28022 143.1347 44.3579 49.5019 33.7255 −6.6773 56.1475 8.2295 54.8241 214.4
0.9236 0.9002 0.9302 0.9866 1.0127 0.9002 0.9628 1.0016 0.9059 0.9895 45549.49 1.2898 18.4275 0.59267 0.30052 310.4204 −29.4561 41.862 58.879 −4.8572 32.9393 8.4604 149.6131 212.5
Fig. 13. Convergence of case 8 (ECHT-DE) for IEEE 30-bus system. Fig. 15. IEEE 57-bus system – voltage profiles of load buses for best solutions of case 11 to case 14.
value of cumulative voltage deviation (VD in Eq. (23)) would be 3 p.u. (i.e. 50 × 0.06 p.u.). Reported VD value in one case is found to be higher than 3 p.u. and the same is marked with a footnote in comparison Table 10. The reference paper pertaining to the case study adopted most common static penalty function approach for constraint handling. Case 11 of minimizing basic objective of fuel cost leads to a value of 41667.82 $/h by SP-DE, the lowest when compared with other valid results of recent studies as depicted in Table 10. Further, network loss is also among the best (14.909 MW) with the recommended settings of control variables by SP-DE algorithm. Study in case 12 of simultaneous minimization of fuel cost and voltage deviation (𝑉 𝐷) by SP-DE achieves fitness value close to the value obtained by DSA (Abaci and Yamacli, 2016) method. Significant improvement in fuel cost is observed (by SFDE) in multi-objective optimization of case 13 where both cost and 𝐿index (max) are minimized. Case 14 is single objective optimization of minimizing voltage deviation. The value reported by SP-DE is much better than the comparable study by APFPA (Mahdad and Srairi, 2016), though with the sacrifice of fuel cost which increases from previous studies following the recommended settings of control variables. Roy and Paul (2015) and Shaheen et al. (2017) have used higher upper limits for shunt VAR compensators (up to 30 MVAr); hence straight comparison with present study is not valid. In general, fitness values
Fig. 14. Convergence of case 9 (SF-DE) for IEEE 30-bus system.
cases except case 13. Similar to system 1 of 30-buses, constraints are all diligently satisfied by the constraint handling methods. Ranges for generator reactive power are derived from Zimmerman et al. (0000) and found to be quite narrow for some generators. The CH techniques duly satisfy the limits as listed values of generator reactive power suggest. Further, the permissible range of load bus voltages for 57-bus system is between [0.94 p.u.−1.06 p.u.]. Theoretically, if all 50-load buses of the system operate at the limits, the maximum possible absolute 94
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Fig. 19. Comparative convergence of case 14 for IEEE 57-bus system.
Fig. 16. Comparative convergence of case 11 for IEEE 57-bus system.
Fig. 17. Convergence of case 12 (SP-DE) for IEEE 57-bus system. Fig. 20. IEEE 118-bus system – voltage profiles of load buses for study cases with SP-DE.
Fig. 18. Convergence of case 13 (SF-DE) for IEEE 57-bus system. Fig. 21. Convergences of case 15 and case 16 (SP-DE) for IEEE 118-bus system.
of the objective functions achieved by SP-DE and SF-DE in all the cases are better than most other reported values in past literatures. Curves in Fig. 15 describe the voltage profiles of load buses for the best solutions of case studies performed for 57-bus system. Clearly all voltages are within allowable limits. However, in a couple of cases many bus voltages are found to be very close to the upper limit. The fact reiterates the need of extra caution in satisfying voltage constraint so that no system bus experiences overvoltage. Convergence curves applying 3-CH techniques are superimposed for case 11 in Fig. 16 which shows reasonably fast convergence to the fitness value. SP-DE appears to have converged faster than the other methods. All the curves are shown starting after a certain number of function
evaluations. It is worthwhile to mention that, during initial phase of the search process, the algorithm (CH technique) looks for feasible solutions. The actual convergence to optimal solution starts when the search process enters the feasible region. For 30-bus system, the CH techniques easily find the feasible region. But possibly due to narrow limits of generator reactive power in 57-bus system, the algorithm takes few hundreds (or sometimes few thousands) of evaluations to find the feasible zone. The convergence curves do not indicate the initial phase of search process when the fitness progression is erratic. Now, for twoobjective cases, convergence curves of only best performing method are drawn in case 12 and case 13 to maintain clarity in diagrams. 95
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 10 Comparison of 3-CH techniques with previous studies for various cases of IEEE 57-bus system. Case #
Algorithm
Fitness
Fuel cost ($/h)
Emission (t/h)
𝑃loss (MW)
VD (p.u.)
L-index (max)
Case 11
ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) DSA (Abaci and Yamacli, 2016) ICBO (Bouchekara et al., 2016b) ARCBBO (Kumar and Premalatha, 2015) APFPA (Mahdad and Srairi, 2016) LTLBO (Ghasemi et al., 2015) MICA-TLA (Ghasemi et al., 2014) DE (Shaheen et al., 2017) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) DSA (Abaci and Yamacli, 2016) MFO (Mohamed et al., 2017) MICA-TLA (Ghasemi et al., 2014) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) DSA (Abaci and Yamacli, 2016) MFO (Mohamed et al., 2017) ECHT-DE SF-DE SP-DE APFPA (Mahdad and Srairi, 2016) KHA (Roy and Paul, 2015) DE (Shaheen et al., 2017)
41670.56 41667.85 41667.82 41673.72 41686.82 41697.33 41686 41628.75a 41679.55 41675.05 41682 41776.48 41775.09 41774.75 41782.80 41775.60 41786.66 42013.08 41699.25 41695.55 41696.54 41703.48 41785.05 41707.66 0.60416 0.59584 0.59267 0.8909 0.5810b 0.5839b
41670.56 41667.85 41667.82 41673.72 41686.82 41697.33 41686 41628.75a 41679.55 41675.05 41682 41694.82 41697.52 41697.50 41714.98 41699.40 41718.87 41959.18 41671.09 41667.53 41668.45 41675.99 41761.22 41680.19 46813.22 45246.02 45549.49 43485.93 42006.44 –
1.36230 1.35816 1.350 1.9526 – – – – – – – 1.3597 1.35769 1.3550 1.9551 – 2.0149 – 1.3609 1.3576 1.35524 1.9188 – 1.9192 1.3379 1.23453 1.2898 – – –
14.9479 14.8864 14.9090 15.0526 – 15.5470 15.3769 14.0470 15.1589 15.0149 – 15.5806 15.5616 15.5897 15.9214 – 16.2189 19.909 15.0275 14.8963 15.012 15.0026 – 15.1026 19.0821 18.4697 18.4275 12.1513 – –
1.50319 1.64209 1.54367 1.5508 1.0833 1.3173 – 3.5571a – 1.6161 – 0.81659 0.77572 0.77253 0.6782 0.7620 0.6780 0.5390 1.56188 1.61174 1.60803 1.7236 1.0573 1.7245 0.60416 0.59584 0.59267 0.8909 0.5810b 0.5839b
0.28886 0.27971 0.28123 0.28392 0.24353 0.27760 – – – – – 0.29198 0.29262 0.29228 0.29533 0.2471 0.29525 – 0.28152 0.28022 0.28092 0.27481 0.2383 0.27467 0.3008 0.30135 0.30052 – 0.2985 –
Case 12
Case 13
Case 14
a
Infeasible solution, constraint on load bus voltage is violated.
b
Higher limits for the shunt compensators have been considered.
As expected, unevenness in individual objectives, especially voltage deviation, along the optimization path is observed for multi-objective cases represented in Figs. 17 and 18. Single objective of minimizing voltage deviation requires higher number of function evaluations for all methods to converge as seen in Fig. 19. It might be because of highly non-linear relationship of bus voltages with numerous control variables of 57-bus system.
5.5. Additional study case, AC1 – handling discrete variables This section of the study considers cost minimization of IEEE-30 bus system considering discrete variables in the optimization process. In DE, discrete variables can easily be accounted for by ‘rounding off ‘operation. Here, both (shunt) capacitor banks and transformer taps are considered discrete variables. The capacitor can be switched in at discrete steps of 200 kVAr (i.e. 0.2 MVAr). Transformer tap is considered in p.u. within range 0.90 p.u. to 1.10 p.u. The possible transformer taps considered in this case study are [0.90–0.92–0.94-. . . −1.08–1.10] p.u. During the search process, the control variables on shunt capacitor ratings and transformer taps are simply rounded off (up to required decimal points) to the nearest multiples of 0.2 (MVAr) and 0.02 (p.u.), respectively. The optimization case on generation cost is run with all the three algorithms i.e. ECHT-DE, SF-DE and SP-DE. SF-DE algorithm achieves the minimum cost value of 800.4352 $/h, which is slightly higher than the best value obtained in same case with all continuous variables. The complete results of all algorithms with control variables and other calculated parameters are provided in Table A.7 in Appendix.
5.4. System 3: IEEE 118-bus test system To verify effectiveness of the proposed method, the large-scale IEEE 118-bus system is considered for study purpose. In general, performance of SP-DE is found to be superior for higher number of variables in constrained optimization problems. Hence, SP-DE algorithm is applied to the system to minimize fuel cost (case 15) and power loss (case 16) independently. A total of 300,000 function evaluations (with execution time of about 40 min) are performed for a study case with population size (𝑁𝑃 ) of 30. The details of calculation of these objective functions are mentioned in Section 3. The control variables and calculated parameters pertaining to the best fitness values obtained in 5-trial runs of each case are provided in Tables A.5 and A.6 in Appendix. Allowable ranges of generator real power (Zimmerman et al., 0000) are explicitly mentioned in the table. Shunt compensator ratings are selectable within 0 to 25 MVAr. Permissible ranges of generator bus voltages and transformer tap settings are within [0.95, 1.10] and [0.90, 1.10], respectively. The optimal value of fuel cost in case 15 is 135055.7 $/h, and that of real power loss in case 16 is 17.6946 MW. The system has a large number of control variables and constraints, and the number of critical constraints on generator reactive power and load bus accounts up to 118 nos. All the constraints are duly satisfied by the algorithm. Load bus voltage profiles for both the study cases are provided in Fig. 20. Convergence characteristics are shown in Fig. 21. Due to presence of huge number of constraints, the algorithm takes about 25 000 function evaluations to find first set of feasible solutions. The convergence to optimal solutions is observed in about 200 000 function evaluations.
6. Conclusion This paper discusses in detail the application and usefulness of three constraint handling techniques – SF, SP and ensemble of these two with DE being the basic search algorithm for optimal power flow (OPF) problem which consists of non-linear objectives and non-linear constraints. A comparative study among the three techniques points out the difficulty in guaranteeing superiority of a CH technique over others under various scenarios of a real-world problem. Ensemble method helps to an extent in achieving near best solution almost in all cases. However, the ensemble method does not always promise most optimal solution and fast convergence. Nonetheless, the importance of an efficient constraint handling technique is paramount. As substantiated in our study, without proper constraint handling method, and especially with static 96
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table A.1 Branch data used for 30-bus system. Branch no.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
Bus no.
R (p.u.)
From
To
1 1 2 3 2 2 4 5 6 6 6 6 9 9 4 12 12 12 12 14 16
2 3 4 4 5 6 6 7 7 8 9 10 11 10 12 13 14 15 16 15 17
0.0192 0.0452 0.0570 0.0132 0.0472 0.0581 0.0119 0.0460 0.0267 0.0120 0 0 0 0 0 0 0.1231 0.0662 0.0945 0.2210 0.0524
X (p.u.)
B/2 (p.u.)
0.0575 0.1652 0.1737 0.0379 0.1983 0.1763 0.0414 0.1160 0.0820 0.0420 0.2080 0.5560 0.2080 0.1100 0.2560 0.1400 0.2559 0.1304 0.1987 0.1997 0.1923
0.0264 0.0204 0.0184 0.0042 0.0209 0.0187 0.0045 0.0102 0.0085 0.0045 0 0 0 0 0 0 0 0 0 0 0
Rating (MVA)
Branch no.
130 130 65 130 130 65 90 70 130 32 65 32 65 65 65 65 32 32 32 16 16
Bus no.
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
From
To
15 18 19 10 10 10 10 21 15 22 23 24 25 25 28 27 27 29 8 6
18 19 20 20 17 21 22 22 23 24 24 25 26 27 27 29 30 30 28 28
R (p.u.)
X (p.u.)
B/2 (p.u.)
Rating (MVA)
0.1073 0.0639 0.0340 0.0936 0.0324 0.0348 0.0727 0.0116 0.1000 0.1150 0.1320 0.1885 0.2544 0.1093 0 0.2198 0.3202 0.2399 0.0636 0.0169
0.2185 0.1292 0.0680 0.2090 0.0845 0.0749 0.1499 0.0236 0.2020 0.1790 0.2700 0.3292 0.3800 0.2087 0.3960 0.4153 0.6027 0.4533 0.2000 0.0599
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0214 0.0065
16 16 32 32 32 32 32 32 16 16 16 16 16 16 65 16 16 16 32 32
Table A.2 Cost and emission coefficients of generators for IEEE 30-bus system (Chaib et al., 2016; Alsac and Stott, 1974). Generator
Bus
𝑎
𝑏
𝑐
𝑑
𝑒
𝛼
𝛽
𝛾
𝜔
μ
𝐺1 𝐺2 𝐺3 𝐺4 𝐺5 𝐺6
1 2 5 8 11 13
0 0 0 0 0 0
2 1.75 1 3.25 3 3
0.00375 0.0175 0.0625 0.00834 0.025 0.025
18 16 14 12 13 13.5
0.037 0.038 0.04 0.045 0.042 0.041
4.091 2.543 4.258 5.326 4.258 6.131
−5.554 −6.047 −5.094 −3.55 −5.094 −5.555
6.49 5.638 4.586 3.38 4.586 5.151
0.0002 0.0005 0.000001 0.002 0.000001 0.00001
2.857 3.333 8 2 8 6.667
Table A.3 Cost coefficients of generators 1 and 2 for IEEE 30-bus system for multi-fuel sources (Chaib et al., 2016). Generator
Bus
𝑎
𝐺1 𝐺2
1 2
55 40
𝑏
𝑐
𝑃𝐺𝑚𝑖𝑛
𝑃𝐺𝑚𝑎𝑥
𝑎
50 20
140 55
82.5 80
𝑖1
0.7 0.3
0.005 0.01
𝑏
𝑐
𝑖1
𝑃𝐺𝑚𝑖𝑛
𝑃𝐺𝑚𝑎𝑥
140 55
200 80
𝑖2
1.05 0.6
0.0075 0.02
𝑖2
Table A.4 Cost and emission coefficients of generators for IEEE 57-bus system (Chaib et al., 2016; Zimmerman et al., 0000). Generator
Bus
𝑎
𝑏
𝑐
𝑑
𝑒
𝛼
𝛽
𝛾
𝜔
μ
𝐺1 𝐺2 𝐺3 𝐺4 𝐺5 𝐺6 𝐺7
1 2 3 6 8 9 12
0 0 0 0 0 0 0
20 40 20 40 20 40 20
0.0775795 0.01 0.25 0.01 0.0222222 0.01 0.0322581
18 16 13.5 18 14 15 12
0.037 0.038 0.041 0.037 0.04 0.039 0.045
4.091 2.543 6.131 3.491 4.258 2.754 5.326
−5.554 −6.047 −5.555 −5.754 −5.094 −5.847 −3.555
6.49 5.638 5.151 6.39 4.586 5.238 3.38
0.0002 0.0005 0.00001 0.0003 0.000001 0.0004 0.002
0.286 0.333 0.667 0.266 0.8 0.288 0.2
help faster convergence and betterment of solutions especially for largescale power systems.
penalty function approach, limits on network parameters may often unknowingly be violated. Violating physical or security constraints of network components may compromise system safety, lead to excessive losses, malfunction and often failure of the component. So, operating the network within defined limits is the pre-requisite for secure and proper operation. Recommended settings by the CH techniques successfully lead to the desired condition of the network. Furthermore, the proposed algorithms in this study perform noticeably better than many other equivalent optimization methods in finding solutions of OPF. Reduction in hourly operation cost has been established almost in all the cases studied under the scope of this literature. In future, the authors propose integration of advanced variants of DE with CH techniques that might
Acknowledgment This project is funded by the National Research Foundation Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) program. Appendix A See Table A.1–A.7. 97
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table A.5 Simulation results of SP-DE algorithm for Case 15 of IEEE 118-bus system. Control variables
Range (MW)
Values (MW)
Control variables
Range (MW)
Values (MW)
Control variables
Values (p.u.)
Control variables
Values (p.u.)
Control variables
Values
𝑃G1 (MW)
30–100
30.0317
𝑃G65
147.3–491
289.6489
𝑉G1 (MW)
0.9871
𝑉G65
1.0623
𝑄C5 (MVAr)
14.0629
𝑃G4
30–100
30.0143
𝑃G66
147.6–492
289.1504
𝑉G4
1.0153
𝑉G66
1.0593
𝑄C34 (MVAr)
0.8808
𝑃G6
30–100
30.0122
𝑃G70
30–100
30.0116
𝑉G6
1.008
𝑉G69
1.0389
𝑄C37 (MVAr)
0.8844
𝑃G8
30–100
30.0052
𝑃G72
30–100
30.0013
𝑉G8
1.0388
𝑉G70
1.0195
𝑄C44 (MVAr)
5.768
𝑃G10
165–550
317.0191
𝑃G73
30–100
30.0033
𝑉G10
1.0494
𝑉G72
1.0191
𝑄C45 (MVAr)
21.5888
𝑃G12
55.5–185
66.8229
𝑃G74
30–100
30.0088
𝑉G12
1.002
𝑉G73
1.0234
𝑄C46 (MVAr)
10.9322
𝑃G15
30–100
30.0162
𝑃G76
30–100
30.0074
𝑉G15
0.9988
𝑉G74
1.0058
𝑄C48 (MVAr)
4.6786
𝑃G18
30–100
30.0063
𝑃G77
30–100
30.0141
𝑉G18
0.9989
𝑉G76
0.9868
𝑄C74 (MVAr)
24.2029
𝑃G19
30–100
30.0495
𝑃G80
173.1–577
350.9989
𝑉G19
0.9986
𝑉G77
1.013
𝑄C79 (MVAr)
23.8787
𝑃G24
30–100
30.0111
𝑃G85
30–100
30.0087
𝑉G24
1.015
𝑉G80
1.0218
𝑄C82 (MVAr)
23.5807
𝑃G25
96–320
152.1726
𝑃G87
31.2–104
31.2015
𝑉G25
1.0298
𝑉G85
1.0242
𝑄C83 (MVAr)
20.4897
𝑃G26
124.2–414
220.8106
𝑃G89
212.1–707
379.9452
𝑉G26
1.0744
𝑉G87
1.0432
𝑄C105 (MVAr)
13.4731
𝑃G27
30–100
30.0364
𝑃G90
30–100
30.0443
𝑉G27
1.0045
𝑉G89
1.0274
𝑄C107 (MVAr)
1.936
𝑃G31
32.1–107
32.1004
𝑃G91
30–100
30.021
𝑉G31
0.9992
𝑉G90
1.0062
𝑄C110 (MVAr)
18.1676
𝑃G32
30–100
30.0143
𝑃G92
30–100
30.0162
𝑉G32
1.0045
𝑉G91
1.0074
𝑇8 (p.u.)
1.0148
𝑃G34
30–100
30.0024
𝑃G99
30–100
30.0027
𝑉G34
1.0141
𝑉G92
1.0153
𝑇32 (p.u.)
1.0978
𝑃G36
30–100
30.0132
𝑃G100
105.6–352
177.1013
𝑉G36
1.0105
𝑉G99
1.0182
𝑇36 (p.u.)
1.0348
𝑃G40
30–100
30.0216
𝑃G103
42–140
42.0053
𝑉G40
1.0001
𝑉G100
1.0187
𝑇51 (p.u.)
1.0107
𝑃G42
30–100
30.038
𝑃G104
30–100
30.0088
𝑉G42
1.0081
𝑉G103
1.0146
𝑇93 (p.u.)
0.9946
𝑃G46
35.7–119
35.7003
𝑃G105
30–100
30.0022
𝑉G46
1.0316
𝑉G104
1.0067
𝑇95 (p.u.)
1.0095
𝑃G49
91.2–304
162.3848
𝑃G107
30–100
30.013
𝑉G49
1.0429
𝑉G105
1.0063
𝑇102 (p.u.)
0.9771
𝑃G54
44.4–148
44.6599
𝑃G110
30–100
30.0043
𝑉G54
1.0217
𝑉G107
0.9993
𝑇107 (p.u.)
0.9715
𝑃G55
30–100
30.0485
𝑃G111
40.8–136
40.8014
𝑉G55
1.0215
𝑉G110
1.0147
𝑇127 (p.u.)
1.0214
𝑃G56
30–100
30.0079
𝑃G112
30–100
30.0166
𝑉G56
1.0215
𝑉G111
1.0247
Fuel cost ($/h)
135055.7
𝑃G59
76.5–255
125.3306
𝑃G113
30–100
30.0223
𝑉G59
1.0424
𝑉G112
1.0046
𝑃loss (MW)
60.9596
𝑃G61
78–260
124.1197
𝑃G116
30–100
30.0052
𝑉G61
1.0496
𝑉G113
1.0057
VD (p.u.)
1.0715
𝑃G62
30–100
30.0168
𝑉G62
1.0462
𝑉G116
1.0592
𝑃G69 (MW) [0 – 805.2]
370.4283
Values (MW)
Control variables
Values (p.u.)
Control variables
Values (p.u.)
Control variables
Values
Table A.6 Simulation results of SP-DE algorithm for Case 16 of IEEE 118-bus system. Control variables
Range (MW)
Values (MW)
Control variables
Range (MW)
𝑃G1 (MW)
30–100
67.1484
𝑃G65
147.3–491
147.3227
𝑉G1 (MW)
0.9888
𝑉G65
1.0134
𝑄C5 (MVAr)
9.3906
𝑃G4
30–100
30.0136
𝑃G66
147.6–492
147.606
𝑉G4
1.005
𝑉G66
1.047
𝑄C34 (MVAr)
0.322
𝑃G6
30–100
30.95
𝑃G70
30–100
30.1083
𝑉G6
1.0009
𝑉G69
1.0177
𝑄C37 (MVAr)
0.0376
𝑃G8
30–100
30.0116
𝑃G72
30–100
30.0151
𝑉G8
1.029
𝑉G70
1.0149
𝑄C44 (MVAr)
3.1065
𝑃G10
165–550
165.0121
𝑃G73
30–100
30.029
𝑉G10
1.04
𝑉G72
1.0215
𝑄C45 (MVAr)
21.2411
𝑃G12
55.5–185
136.4669
𝑃G74
30–100
95.4903
𝑉G12
1
𝑉G73
1.02
𝑄C46 (MVAr)
5.8559
𝑃G15
30–100
84.631
𝑃G76
30–100
99.8349
𝑉G15
1.0008
𝑉G74
1.0127
𝑄C48 (MVAr)
8.9501
𝑃G18
30–100
32.5503
𝑃G77
30–100
99.9777
𝑉G18
0.9997
𝑉G76
0.9992
𝑄C74 (MVAr)
23.8294
𝑃G19
30–100
58.8217
𝑃G80
173.1–577
289.4337
𝑉G19
0.9992
𝑉G77
1.0079
𝑄C79 (MVAr)
24.9914
𝑃G24
30–100
30.0197
𝑃G85
30–100
31.0829
𝑉G24
1.0164
𝑉G80
1.0131
𝑄C82 (MVAr)
24.8441
𝑃G25
96–320
96.0079
𝑃G87
31.2–104
31.2094
𝑉G25
1.0239
𝑉G85
1.0124
𝑄C83 (MVAr)
14.0215
𝑃G26
124.2–414
124.2228
𝑃G89
212.1–707
212.105
𝑉G26
1.0443
𝑉G87
1.0332
𝑄C105 (MVAr)
11.5006
𝑃G27
30–100
49.1688
𝑃G90
30–100
99.7784
𝑉G27
1.0062
𝑉G89
1.0181
𝑄C107 (MVAr)
20.4597
𝑃G31
32.1–107
61.5931
𝑃G91
30–100
30.0022
𝑉G31
1.0048
𝑉G90
1.0096
𝑄C110 (MVAr)
17.3712
𝑃G32
30–100
39.6414
𝑃G92
30–100
30.1177
𝑉G32
1.0058
𝑉G91
1.0128
𝑇8 (p.u.)
1.0224
𝑃G34
30–100
63.2896
𝑃G99
30–100
38.6397
𝑉G34
1.0055
𝑉G92
1.0093
𝑇32 (p.u.)
1.0727
𝑃G36
30–100
53.3235
𝑃G100
105.6–352
105.6473
𝑉G36
1.0015
𝑉G99
1.0114
𝑇36 (p.u.)
1.0129
𝑃G40
30–100
99.9392
𝑃G103
42–140
42.1357
𝑉G40
1.0087
𝑉G100
1.0116
𝑇51 (p.u.)
1.0044
𝑃G42
30–100
99.9191
𝑃G104
30–100
30.8949
𝑉G42
1.0175
𝑉G103
1.0138
𝑇93 (p.u.)
0.9599
𝑃G46
35.7–119
82.1758
𝑃G105
30–100
52.9999
𝑉G46
1.0371
𝑉G104
1.0117
𝑇95 (p.u.)
0.9583
𝑃G49
91.2–304
141.0635
𝑃G107
30–100
57.6295
𝑉G49
1.0395
𝑉G105
1.0123
𝑇102 (p.u.)
0.9431
𝑃G54
44.4–148
147.8306
𝑃G110
30–100
30.1049
𝑉G54
1.0519
𝑉G107
1.012
𝑇107 (p.u.)
0.9677 0.9847
𝑃G55
30–100
79.5092
𝑃G111
40.8–136
40.8554
𝑉G55
1.052
𝑉G110
1.017
𝑇127 (p.u.)
𝑃G56
30–100
99.8073
𝑃G112
30–100
51.6249
𝑉G56
1.0509
𝑉G111
1.0247
Fuel cost ($/h)
155724.9
𝑃G59
76.5–255
251.0927
𝑃G113
30–100
30.0672
𝑉G59
1.05
𝑉G112
1.014
𝑃loss (MW)
17.6946
𝑃G116
30–100
75.151
𝑃G61
78–260
78.2118
𝑃G62
30–100
64.477
𝑉G61
1.0493
𝑉G113
1.0075
VD (p.u.)
0.8663
𝑉G62
1.0473
𝑉G116
1.0101
𝑃G69 (MW) [0 – 805.2]
2.9323
98
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100 Table A.7 Simulation results of study case AC1: handling discrete variables. Parameters
Range
ECHT-DE
SF-DE
SP-DE
𝑃G2 (MW) 𝑃G5 (MW) 𝑃G8 (MW) 𝑃G11 (MW) 𝑃G13 (MW) 𝑉1 (p.u.) 𝑉2 (p.u.) 𝑉5 (p.u.) 𝑉8 (p.u.) 𝑉11 (p.u.) 𝑉13 (p.u.) 𝑄c10 (MVAr) 𝑄c12 (MVAr) 𝑄c15 (MVAr) 𝑄c17 (MVAr) 𝑄c20 (MVAr) 𝑄c21 (MVAr) 𝑄c23 (MVAr) 𝑄c24 (MVAr) 𝑄c29 (MVAr) 𝑇11 (p.u.) 𝑇12 (p.u.) 𝑇15 (p.u.) 𝑇36 (p.u.) Fuel cost ($/h) Emission (t/h) 𝑃loss (MW) VD (p.u.) L-index (max) 𝑃G1 (MW)
20 – 80 15 – 50 10 – 35 10 – 30 12 – 40 0.95 – 1.10 0.95 – 1.10 0.95 – 1.10 0.95 – 1.10 0.95 – 1.10 0.95 – 1.10 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.90 – 1.10 0.90 – 1.10 0.90 – 1.10 0.90 – 1.10
48.4217 21.5846 21.3143 11.9471 12.0404 1.0821 1.0631 1.032 1.0357 1.0676 1.0533 3.2 1.2 4.4 5 3.4 5 3 5 2.8 1.04 0.92 0.98 0.98 800.4430 0.36601 8.9896 0.8689 0.1392 177.0814
48.5786 21.4729 21.0862 11.91 12.0122 1.0817 1.0627 1.0314 1.0363 1.0767 1.0532 1 1.4 4.4 5 3.6 5 2.8 5 2.8 1.04 0.92 0.98 0.98 800.4352 0.36684 9.0179 0.87519 0.13906 177.3581
48.36621 21.34167 22.08471 11.80553 12.03479 1.083016 1.064776 1.0343 1.039547 1.08065 1.037134 4.2 3.2 5 4.4 3.2 4.6 3.2 4.8 3 1.06 0.92 0.96 0.98 800.4478 0.36510 8.9648 0.88987 0.13860 176.7319
50 – 200
References
Deb, K., 2000. An efficient constraint handling method for genetic algorithms. Comput. Methods Appl. Mech. Eng. 186 (2), 311–338. Derrac, J., García, S., Molina, D., Herrera, F., 2011. A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms. Swarm Evol. Comput. 1 (1), 3–18. Farmani, R., Wright, J.A., 2003. Self-adaptive fitness formulation for constrained optimization. IEEE Trans. Evol. Comput. 7 (5), 445–455. Ghasemi, M., Ghavidel, S., Gitizadeh, M., Akbari, E., 2015. An improved teaching– learning-based optimization algorithm using Lévy mutation strategy for non-smooth optimal power flow. Int. J. Electr. Power Energy Syst. 65, 375–384. Ghasemi, M., Ghavidel, S., Rahmani, S., Roosta, A., Falah, H., 2014. A novel hybrid algorithm of imperialist competitive algorithm and teaching learning algorithm for optimal power flow problem with non-smooth cost functions. Eng. Appl. Artif. Intell. 29, 54–69. Kessel, P., Glavitsch, H., 1986. Estimating the voltage stability of a power system. IEEE Trans. Power Deliv. 1 (3), 346–354. Kumar, A.R., Premalatha, L., 2015. Optimal power flow for a deregulated power system using adaptive real coded biogeography-based optimization. Int. J. Electr. Power Energy Syst. 73, 393–399. Li, C., Yang, S., Nguyen, T.T., Yu, E.L., Yao, X., Jin, Y., Beyer, H.G., Suganthan, P.N., (2008). Benchmark generator for CEC 2009 competition on dynamic optimization. Liang, J.J., Qu, B.Y., Suganthan, P.N., (2013). Problem definitions and evaluation criteria for the CEC 2014 special session and competition on single objective real-parameter numerical optimization. Computational Intelligence Laboratory, Zhengzhou University, Zhengzhou China and Technical Report, Nanyang Technological University, Singapore. Liang, J.J., Runarsson, T.P., Mezura-Montes, E., Clerc, M., Suganthan, P.N., Coello, C.C., Deb, K., 2006. Problem definitions and evaluation criteria for the CEC 2006 special session on constrained real-parameter optimization. J. Appl. Mech. 41 (8). Mahdad, B., Srairi, K., 2016. Security constrained optimal power flow solution using new adaptive partitioning flower pollination algorithm. Appl. Soft Comput. 46, 501–522. Mallick, S., Acharjee, P., Ghoshal, S.P., Thakur, S.S., 2013. Determination of maximum load margin using fuzzy logic. Int. J. Electr. Power Energy Syst. 52, 231–246. Mallipeddi, R., Jeyadevi, S., Suganthan, P.N., Baskar, S., 2012. Efficient constraint handling for optimal reactive power dispatch problems. Swarm Evol. Comput. 5, 28– 36. Mallipeddi, R., Suganthan, P.N., 2010. Ensemble of constraint handling techniques. IEEE Trans. Evol. Comput. 14 (4), 561–579. Mohamed, A.A.A., Mohamed, Y.S., El-Gaafary, A.A., Hemeida, A.M., 2017. Optimal power flow using moth swarm algorithm. Electr. Power Syst. Res. 142, 190–206. Naderi, E., Narimani, H., Fathi, M., Narimani, M.R., 2017. A novel fuzzy adaptive configuration of particle swarm optimization to solve large-scale optimal reactive power dispatch. Appl. Soft Comput. 53, 441–456.
Abaci, K., Yamacli, V., 2016. Differential search algorithm for solving multi-objective optimal power flow problem. Int. J. Electr. Power Energy Syst. 79, 1–10. Acharjee, P., 2012. Identification of maximum loadability limit and weak buses using security constraint genetic algorithm. Int. J. Electr. Power Energy Syst. 36 (1), 40–50. Acharjee, P., Mallick, S., Thakur, S.S., Ghoshal, S.P., 2011. Detection of maximum loadability limits and weak buses using Chaotic PSO considering security constraints. Chaos Solut. Fractals 44 (8), 600–612. Alsac, O., Stott, B., 1974. Optimal load flow with steady-state security. IEEE Trans. Power Apparatus Syst. (3), 745–751. Ayan, K., Kılıç, U., Baraklı, B., 2015. Chaotic artificial bee colony algorithm based solution of security and transient stability constrained optimal power flow. Int. J. Electr. Power Energy Syst. 64, 136–147. Bhattacharya, A., Chattopadhyay, P.K., 2011. Application of biogeography-based optimisation to solve different optimal power flow problems. IET Gener. Transm. Distribution 5 (1), 70–80. Bhowmik, A.R., Chakraborty, A.K., 2015. Solution of optimal power flow using nondominated sorting multi objective opposition based gravitational search algorithm. Int. J. Electr. Power Energy Syst. 64, 1237–1250. Biswas, P.P., Suganthan, P.N., Amaratunga, G.A., 2017. Optimal power flow solutions incorporating stochastic wind and solar power. Energy Convers. Manag. 148, 1194– 1207. Bouchekara, H.R.E.H., Chaib, A.E., Abido, M.A., 2016a. Multiobjective optimal power flow using a fuzzy based grenade explosion method. Energy Syst. 7 (4), 699–721. Bouchekara, H.R.E.H., Chaib, A.E., Abido, M.A., El-Sehiemy, R.A., 2016b. Optimal power flow using an improved colliding bodies optimization algorithm. Appl. Soft Comput. 42, 119–131. Cai, H.R., Chung, C.Y., Wong, K.P., 2008. Application of differential evolution algorithm for transient stability constrained optimal power flow. IEEE Trans. Power Syst. 23 (2), 719–728. Chaib, A.E., Bouchekara, H.R.E.H., Mehasni, R., Abido, M.A., 2016. Optimal power flow with emission and non-smooth cost functions using backtracking search optimization algorithm. Int. J. Electr. Powerand Energy Syst. 81, 64–77. Chen, Q., Liu, B., Zhang, Q., Liang, J.J., Suganthan, P.N., Qu, B.Y., (2014). Problem definition and evaluation criteria for CEC 2015 special session and competition on bound constrained single-objective computationally expensive numerical optimization. Computational Intelligence Laboratory, Zhengzhou University, China and Nanyang Technological University, Singapore, Tech. Rep. Coello, C.A.C., 2002. Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: A survey of the state of the art. Comput. Methods Appl. Mech. Eng. 191 (11), 1245–1287. Daryani, N., Hagh, M.T., Teimourzadeh, S., 2016. Adaptive group search optimization algorithm for multi-objective optimal power flow problem. Appl. Soft Comput. 38, 1012–1024. 99
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Niknam, T., Azizipanah-Abarghooee, R., Narimani, M.R., 2012. A new multi objective optimization approach based on TLBO for location of automatic voltage regulators in distribution systems. Eng. Appl. Artif. Intell. 25 (8), 1577–1588. Pandiarajan, K., Babulal, C.K., 2016. Fuzzy harmony search algorithm based optimal power flow for power system security enhancement. Int. J. Electr. Power Energy Syst. 78, 72–79. Pulluri, H., Naresh, R., Sharma, V., 2016. A solution network based on stud krill herd algorithm for optimal power flow problems. Soft. Comput. 1–18. Pulluri, H., Naresh, R., Sharma, V., 2017. Application of stud krill herd algorithm for solution of optimal power flow problems. Int. Trans. Electr. Energy Syst. 27 (6). Qin, A.K., Huang, V.L., Suganthan, P.N., 2009. Differential evolution algorithm with strategy adaptation for global numerical optimization. IEEE Trans. Evol. Comput. 13 (2), 398–417. Rabiee, A., Parniani, M., 2013. Voltage security constrained multi-period optimal reactive power flow using benders and optimality condition decompositions. IEEE Trans. Power Syst. 28 (2), 696–708. Reddy, S.S., Rathnam, C.S., 2016. Optimal power flow using glowworm swarm optimization. Int. J. Electr. Power Energy Syst. 80, 128–139. Roy, P.K., Paul, C., 2015. Optimal power flow using krill herd algorithm. Int. Trans. Electr. Energy Syst. 25 (8), 1397–1419.
Roy, R., Jadhav, H.T., 2015. Optimal power flow solution of power system incorporating stochastic wind power using Gbest guided artificial bee colony algorithm. Int. J. Electr. Power Energy Syst. 64, 562–578. Sedighizadeh, M., Sarvi, M., Naderi, E., 2011. Multi objective optimal power flow with FACTS devices using shuffled frog leaping algorithm. Int. Rev. Electr. Eng. 6 (4). Shaheen, A.M., El-Sehiemy, R.A., Farrag, S.M., 2016. Solving multi-objective optimal power flow problem via forced initialised differential evolution algorithm. IET Gener. Transm. Distribution 10 (7), 1634–1647. Shaheen, A.M., Farrag, S.M., El-Sehiemy, R.A., 2017. MOPF solution methodology. IET Gener. Transm. Distribution 11 (2), 570–581. Storn, R., Price, K., 1997. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optim. 11 (4), 341–359. Tessema, B., Yen, G.G., 2006. A self-adaptive penalty function based algorithm for constrained optimization. In: IEEE Congress on Evolutionary Computation, 2006. CEC 2006. IEEE, pp. 246–253. Wolpert, D.H., Macready, W.G., 1997. No free lunch theorems for optimization. IEEE Trans. Evol. Comput. 1 (1), 67–82. Zimmerman, R.D., Murillo-Sánchez, C.E., Thomas, R.J., Matpower. Available at: http: //www.pserc.cornell.edu/matpower.
100
Contents lists available at ScienceDirect
Engineering Applications of Artificial Intelligence journal homepage: www.elsevier.com/locate/engappai
Optimal power flow solutions using differential evolution algorithm integrated with effective constraint handling techniques Partha P. Biswas a, *, P.N. Suganthan a , R. Mallipeddi b, *, Gehan A.J. Amaratunga c a b c
School of Electrical and Electronic Engineering, Nanyang Technological University, Singapore Kyungpook National University, Daegu, South Korea Department of Engineering, University of Cambridge, UK
a r t i c l e
i n f o
Keywords: Optimal power flow Power loss Emission Differential evolution (DE) Constraint handling (CH) Ensemble of constraint handling techniques (ECHT)
a b s t r a c t Optimal power flow (OPF) is a highly non-linear complex optimization problem where the steady state parameters of an electrical network need to be determined for its economical and efficient operation. The complexity of the problem escalates with ubiquitous presence of constraints in the problem. Solving OPF remains a popular but challenging task among power system researchers. In last couple of decades, numerous evolutionary algorithms (EAs) have been applied to find optimal solutions with different objectives of OPF. However, the search method adopted by EAs is unconstrained. An extensively used methodology to discard infeasible solutions found during the search process is the static penalty function approach. The process requires appropriate selection of penalty coefficients decided largely by tedious trial and error method. This paper presents performance evaluation of proper constraint handling (CH) techniques — superiority of feasibly solutions (SF), self-adaptive penalty (SP) and an ensemble of these two constraint handling techniques (ECHT) with differential evolution (DE) being the basic search algorithm, on the problem of OPF. The methods are tested on standard IEEE 30, IEEE 57 and IEEE 118-bus systems for several OPF objectives such as cost, emission, power loss, voltage stability etc. Single objective and weighted sum multi-objective cases of OPF are studied under the scope of this literature. Simulation results are analyzed and compared with most recent studies on the problem. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction Since its inception more than half a century ago, the optimal power flow (OPF) remains a widely cultivated topic among power system research communities across the globe due to the intriguing multifaceted challenges it poses. The OPF is formulated as single or multiobjective problem of minimizing fuel cost, emission, transmission loss, voltage deviation etc. with constraints to be satisfied on generator capability, line capacity, bus voltage and power flow balance. Output of OPF program helps to determine the optimal operating state of a power system and the corresponding settings of control variables for economic and secure operation (Cai et al., 2008). Major control variables refer to the generated real power and generator bus voltages of the network. The latter controls the reactive power flow which is further compensated by adding capacitor banks of appropriate ratings to the network feeding usually the inductive loads. Voltage vectors of load buses and complex power flows in the lines determined during the process of optimization
represent the system optimal operating state that would result in single or multiple objectives of the network being largely fulfilled. In summary, OPF involves intricate calculations with multiple variables and finding optimal solution while satisfying all constraints simultaneously is the most difficult part one encounters. In earlier days, in use were classical numerical optimization methods which suffered from convexity, assumption of continuity and normally employed a gradient based search that converged to local optima. Revolution in numerical optimization introduced several evolutionary algorithms and techniques in last few decades. Most of these methods can successfully overcome the problem of premature convergence and are able to explore the search space in pursuit of global optima. OPF problem has seen application of numerous such evolutionary algorithms. A few standard objectives of OPF were optimized in Abaci and Yamacli (2016) for IEEE bus systems using differential search algorithm (DSA), an effective algorithm for real-valued numerical optimization problems. Daryani et al. (2016) improved standard group
* Corresponding authors.
E-mail addresses: [email protected] (P.P. Biswas), [email protected] (P.N. Suganthan), [email protected] (R. Mallipeddi), [email protected] (G.A.J. Amaratunga). https://doi.org/10.1016/j.engappai.2017.10.019 Received 24 June 2017; Received in revised form 5 October 2017; Accepted 25 October 2017 0952-1976/© 2017 Elsevier Ltd. All rights reserved.
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
members would result in restricted scope of exploration and exploitation by the algorithm during the search process. System security constraints gained importance in determining maximum loading limits of selected networks in some literatures using genetic algorithm (GA) (Acharjee, 2012), fuzzy logic (Mallick et al., 2013) and chaotic PSO (Acharjee et al., 2011). In optimal reactive power flow problem, voltage security constraint is considered in Rabiee and Parniani (2013). However, state-of-the-art constraint handling (CH) techniques in OPF problem remain largely untested. In present study, superiority of feasible solutions (SF), self-adaptive penalty (SP) methods and an ensemble of these two constraint handling techniques (ECHT) are applied to OPF problem with differential evolution (DE) (Storn and Price, 1997) as the basic search algorithm. Performance of these CH techniques is assessed, results are statistically compared and analyzed in detail. It is worthwhile to note that OPF problem has both equality and inequality constraints. Equality constraints are the power balance equations where both active and reactive power generated in the network must be equal to the demand and losses in the network. Convergence of power flow to a solution ensures that the power balance equations are automatically satisfied. Inequality constraints on slack generator power, reactive power output of the generators, bus voltage limits and line capacities need special attention. SF and SP methods of constraint handling work differently on their own way in handling inequality constraints. The motivation behind ECHT is to assess performance of both the constraint handling techniques by the combinatorial algorithm (i.e. ECHT) as one technique might work better than the other on some problems but not on all problems. Instead of user finding the superior performer out of the two for a particular problem, ECHT takes the burden to apply both the techniques and output close to, if not the best results obtained by one of the constituting techniques. The constraint handling (CH) methods are all successfully implemented in OPF problem on standard IEEE 30, IEEE 57 and IEEE 118-bus systems in our research presented in this paper. Generation cost, emission, real power loss and voltage stability of the network are all individually (as single objective) optimized. Study cases with multi-fuel options for the generators and valve-point loading effect are also considered. Although the study cases mostly consider all continuous variables, handling of discrete variables is also discussed and a study-case results with discrete variables are appended to the result section. Weighted sum approach is adopted in multi-objective optimization for selected study cases of both 30-bus and 57-bus systems. Simulation results are compared and critically investigated against constraint violation with most recent studies on optimal power flow (OPF) found in the literatures. A summary of contributions of this work can be listed pointwise as below:
search optimization algorithm with adaptation to develop adaptive group search optimization (AGSO) to perform study on OPF. Chaib et al. (2016) applied backtracking search optimization algorithm (BSA) to perform OPF calculation with multi-fuel options and valve-point effect in thermal generators. Improved colliding bodies optimization (ICBO) algorithm (Bouchekara et al., 2016b) increased the number of colliding bodies in each iteration to enhance performance of the algorithm when applied to the problem of OPF. Mohamed et al. (2017) applied moth swarm algorithm (MSA) on numerous objectives of OPF for various bus systems to exhibit fast execution time and quick convergence of the algorithm. Chaos theory was incorporated in artificial bee colony to form chaotic artificial bee colony (CABC) in Ayan et al. (2015) and security constrained OPF was solved. Agent based gravitational search algorithm (GSA), where masses of agents play vital roles in guiding the search process, was adopted in Bhowmik and Chakraborty (2015) to solve OPF. To enhance exploration capability and population diversity of biogeography based optimization algorithm (BBO), adaptive real coded BBO (ARCBBO) was suggested in Kumar and Premalatha (2015) and applied to OPF problem. Dynamic adjustment of control parameters in adaptive partitioning flower pollination algorithm (APFPA) (Mahdad and Srairi, 2016) showed higher convergence speed and better accuracy than FPA for OPF solutions. In fuzzy harmony search algorithm (FHSA) (Pandiarajan and Babulal, 2016), effect of fuzzy based automatic adjustment of pitch adjustment rate and bandwidth of harmony search algorithm was studied in applying to the problem of OPF. Fuzzy based adaptive configuration of particle swarm optimization was established in Naderi et al. (2017) to solve a sub-problem of OPF. In solving OPF, Levy mutation strategy for teaching–learning based optimization (LTLBO) (Ghasemi et al., 2015) introduced a Levy mutation operator to enhance exploration at early search stages. Krill herd algorithm (KHA) (Roy and Paul, 2015) and its variant stud krill herd algorithm (SKH) (Pulluri et al., 2017, 2016) have also been popular in finding OPF solutions. Grenade explosion method (GEM) (Bouchekara et al., 2016a), glowworm swarm optimization (GSO) (Reddy and Rathnam, 2016) and Gbest guided artificial bee colony (GABC) (Roy and Jadhav, 2015) have also shown promising results for OPF as reported in respective literatures. Modified imperialist competitive algorithm (MICA) and teaching–learning algorithm (TLA) were hybridized (MICA-TLA) in Ghasemi et al. (2014) to improve local search and convergence of original ICA algorithm in finding OPF solutions. Multi-objective formulation and solution of OPF are found in Shaheen et al. (2016, 2017), Sedighizadeh et al. (2011) and Niknam et al. (2012). Our study focuses on single and weighted-sum multiobjective optimization cases of OPF. As evident from the literature review, either basic or the improved version of an algorithm was tried for real-parameter optimization problem of OPF. Each of the methods has its own strengths and weaknesses, and as depicted in No Free Lunch theorem (Wolpert and Macready, 1997), no single optimization method is capable of solving all types of real-world problems in an effective manner. Differential evolution (DE) and its variants are found to be among the top performing algorithms for single objective real parameter optimization as reflected in the results of CEC competitions (Liang et al., 2006; Li et al., 2008; Liang et al., 2013; Chen et al., 2014). We adopt DE as the basic search algorithm for optimization of OPF objectives. On the aspect constraint handling in OPF, literature (Daryani et al., 2016) defined security objective entwining a couple of inequality constraints as one of the multiple objectives to be minimized. Rest all literatures adopted either static penalty function approach or straightaway discarded the population members that led to constraint violation. As mentioned beforehand, penalty function approach is sensitive to selection of penalty coefficient. Small penalty coefficient over-explores the infeasible region, thus delaying the process of finding feasible solutions, and may prematurely converge to an infeasible solution. On the other hand, large penalty coefficient may not explore the infeasible region properly, thereby resulting in untimely convergence (Mallipeddi et al., 2012). In second approach, removal of infeasible population
∙ Applying CH techniques, superiority of feasible solutions (SF) and self-adaptive penalty (SP) individually to handle constraints while optimizing various objectives of OPF using differential evolution (DE). ∙ Applying ensemble method of CH techniques (ECHT) with SF and SP to handle constraints while optimizing various objectives of OPF using DE. ∙ Statistically compare the results and analyze the performance of SF-DE, SP-DE and ECHT-DE. ∙ Comparison of results of current study with most recent studies on OPF with similar experimental set-up. ∙ Critical analysis of results especially against constraint violation. The organization of rest of the paper is done in following way. Section 2 includes a review of mathematical model including applicable constraints pertaining to OPF in the network. In Section 3, objectives and case studies performed for all the bus test systems are explained with useful numerical values. Description and application of constraint handling (CH) techniques are elaborated in Section 4. Section 5 discusses the simulation results and comparison followed by concluding remarks in Section 6. 82
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
2. Optimal power flow mathematical formulation
2.3.2. Inequality constraints The inequality constraints in OPF reflect the operating limits of the equipment present in the power system and also the limits imposed on lines and load buses to guarantee system security. (a) Generator constraints:
The OPF is a non-linear, non-convex optimization problem and it minimizes certain objectives in the power system subject to many equality and inequality constraints. Mathematically, OPF is represented as:
subject to ∶
𝑉𝐺𝑚𝑖𝑛 ≤ 𝑉𝐺𝑖 ≤ 𝑉𝐺𝑚𝑎𝑥 ∀ 𝑖 ∈ 𝑁𝐺
(6)
𝑃𝐺𝑚𝑖𝑛 ≤ 𝑃𝐺𝑖 ≤ 𝑃𝐺𝑚𝑎𝑥 ∀ 𝑖 ∈ 𝑁𝐺
(7)
𝑚𝑎𝑥 𝑄𝑚𝑖𝑛 ∀ 𝑖 ∈ 𝑁𝐺. 𝐺 ≤ 𝑄𝐺𝑖 ≤ 𝑄𝐺
(8)
𝑖
(1)
Minimize ∶ 𝑓 (𝑥, 𝑢)
𝑖
𝑖
𝑔(𝑥, 𝑢) ≤ 0 ℎ(𝑥, 𝑢) = 0
𝑖
𝑖
𝑖
(b) Transformer constraints:
where, 𝑢 is the vector of control or independent variables, 𝑥 is the vector of state or dependent variables. 𝑓 (𝑥, 𝑢): objective functions of OPF, g (𝑥, 𝑢): set of inequality constraints, ℎ (𝑥, 𝑢): set of equality constraints.
𝑇𝑗𝑚𝑖𝑛 ≤ 𝑇𝑗 ≤ 𝑇𝑗𝑚𝑎𝑥 ∀ 𝑗 ∈ 𝑁𝑇 .
(9)
(c) Shunt compensator constraints: 2.1. Control (independent) variables 𝑚𝑎𝑥 𝑄𝑚𝑖𝑛 ∀ 𝑘 ∈ 𝑁𝐶. 𝐶 ≤ 𝑄 𝐶𝑘 ≤ 𝑄 𝐶 𝑘
The set of variables that can control the power flow in the network is represented in vector form as: 𝑢 = [𝑃𝐺2 … 𝑃𝐺𝑁𝐺 , 𝑉𝐺1 ..𝑉𝐺𝑁𝐺 , 𝑄𝐶1 … 𝑄𝐶𝑁𝐶 , 𝑇1 … 𝑇𝑁𝑇 ]
(2)
(12)
∀ 𝑞 ∈ 𝑛𝑙.
To assess performance of various constraint handling (CH) techniques, several case studies with single and multi-objectives have been performed for standard IEEE 30, 57 and 118-bus test systems. Table 1 includes a summary of major components (viz. generators, transformers, shunt compensators etc.) installed in the test systems and other relevant useful parameters. In 30-bus and 57-bus systems, bus 1 is the swing bus or slack bus, often termed as 𝑉 𝛿 bus. Bus 69 is the swing bus in 118-bus system. In power flow study, the role of swing bus is to balance active and reactive power in the system satisfying power balance Eqs. (4) and (5). For convenience, voltage magnitude (𝑉 ) of swing bus is considered as 1 p.u. while the voltage angle (𝛿) as 0 degree. All other bus voltages and their angles are expressed w.r.t the swing bus, the values of which are obtained as the output of load flow study. Formulation of study cases for the bus test systems is subsequently described.
(3)
3.1. IEEE 30-bus system
2.3. Constraints
A total of ten study cases studies are performed for the system. First six case studies minimize single objective function of OPF. The remaining cases are for multi-objective optimizations which are converted to single objectives with weight factors as in many previous studies and reproduced herein.
As mentioned before, OPF problem has both equality and inequality constraints which are to be satisfied. The constraints are segregated and provided herein. 2.3.1. Equality constraints In OPF, power balance equations are the equality constraints and those are represented as: ( ) 𝑉𝑗 [𝐺𝑖𝑗 cos 𝛿𝑖𝑗 + 𝐵𝑖𝑗 sin(𝛿𝑖𝑗 )] = 0 ∀ 𝑖 ∈ 𝑁𝐵
(4)
( ) 𝑉𝑗 [𝐺𝑖𝑗 sin 𝛿𝑖𝑗 − 𝐵𝑖𝑗 cos(𝛿𝑖𝑗 )] = 0 ∀ 𝑖 ∈ 𝑁𝐵
(5)
3.1.1. Case 1: Minimization of fuel cost This is the most basic objective function of OPF studied almost in all literatures. The association between fuel cost ($/h) and generated power (MW) is approximately given by the quadratic relationship and hence the objective function to be minimized is described as:
𝑗=1 𝑁𝐵 ∑
(11)
3. Objective functions and study cases
where, 𝑃𝐺1 is the generator active power at slack (or swing) bus, 𝑄𝐺𝑖 is the reactive power of generator connected to bus 𝑖, 𝑉𝐿𝑝 is the bus voltage of 𝑝th load bus (PQ bus) and line loading of 𝑞th line is given by 𝑆𝑙𝑞 . 𝑁𝐿 and 𝑛𝑙 are the number of load buses and transmission lines respectively.
𝑄𝐺𝑖 − 𝑄𝐷𝑖 − 𝑉𝑖
𝑝
𝑆𝑙𝑚𝑎𝑥 𝑞
The control variables among the inequality constraints are self-limiting. The optimization algorithm selects a feasible value for each such variable within the defined range. Effective techniques for handling of inequality constraints pertaining to state or dependent variables are discussed in Section 4.
The state of power system is defined by the state variables which can be expressed by vector 𝑥 as:
𝑁𝐵 ∑
≤ 𝑉𝐿𝑝 ≤ 𝑉𝐿𝑚𝑎𝑥 ∀ 𝑝 ∈ 𝑁𝐿
𝑆𝑙𝑞 ≤
2.2. State (dependent) variables
𝑃𝐺𝑖 − 𝑃𝐷𝑖 − 𝑉𝑖
(10)
(d) Security constraints: 𝑉𝐿𝑚𝑖𝑛 𝑝
where, 𝑃𝐺𝑖 is the 𝑖th bus generator active power (except swing generator). Selection of bus 1 as swing bus is representative only and swing bus can be any one of the generator buses. 𝑉𝐺𝑖 is the voltage magnitude at 𝑖th PV bus (generator bus), 𝑇𝑗 is the 𝑗th branch transformer tap, 𝑄𝐶𝑘 is the shunt compensation at 𝑘th bus. 𝑁𝐺, 𝑁𝐶 and 𝑁𝑇 are the number of generators, shunt VAR compensators and transformers respectively. A control variable can assume any value within its range. In reality, transformer taps are not continuous. However, the tap settings expressed here are in p.u. and absolute value of voltage is not accounted for. Hence for study purpose and to compare with past reported results, all control variables including tap settings are considered continuous for most of the study cases. Discrete steps for transformers and shunt capacitors are accounted only in one special study case.
𝑥 = [𝑃𝐺1 , 𝑉𝐿1 … .𝑉𝐿𝑁𝐿 , 𝑄𝐺1 … 𝑄𝐺𝑁𝐺 , 𝑆𝑙1 … 𝑆𝑙𝑛𝑙 ]
𝑘
𝑗=1
𝑓 (𝑥, 𝑢) =
𝑁𝐺 ∑ 𝑖=1
where 𝛿𝑖𝑗 = 𝛿𝑖 − 𝛿𝑗 , is the difference in voltage angles between bus 𝑖 and bus 𝑗, 𝑁𝐵 is the number of buses, 𝑃𝐷 and 𝑄𝐷 are active and reactive load demands, respectively. 𝐺𝑖𝑗 is the transfer conductance and 𝐵𝑖𝑗 is the susceptance between bus 𝑖 and bus 𝑗, respectively.
𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2
𝑖
(13)
where 𝑎𝑖 , 𝑏𝑖 , 𝑐𝑖 are the cost coefficients of the 𝑖th generator producing output power 𝑃𝐺𝑖 . Cost coefficients of IEEE 30-bus system generators are listed in Table A.2 in Appendix. 83
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 1 Summary of systems under study. Items
Buses Branches Generators
IEEE 30-bus system
IEEE 57-bus system
IEEE 118-bus system
Quantity
Details
Quantity
Details
Quantity
Details
30 41 6
Alsac and Stott (1974) Table A.1 (Appendix) Buses: 1 (swing), 2, 5, 8, 11 and 13
57 80 7
Zimmerman et al. (0000) Zimmerman et al. (0000) Buses: 1 (swing), 2, 3, 6, 8, 9 and 12
118 186 54
Zimmerman et al. (0000) Zimmerman et al. (0000) Buses: 1, 4, 6, 8, 10, 12, 15, 18, 19, 24, 25, 26, 27, 31, 32, 34, 36, 40, 42, 46, 49, 54, 55, 56, 59, 61, 62, 65, 66, 69 (swing), 70, 72, 73, 74, 76, 77, 80, 85, 87, 89, 90, 91, 92, 99, 100, 103, 104, 105, 107, 110, 111, 112, 113 and 116 Buses: 5, 34, 37, 44, 45, 46, 48, 74, 79, 82, 83, 105, 107 and 110
Shunt VAR compensation
9
Transformer with tap changer
4
Control variables Connected load Load bus voltage range allowed
24 – 24
Buses: 10, 12, 15, 17, 20, 21, 23, 24 and 29 Branches: 11, 12, 15 and 36
– 283.4 MW, 126.2 MVAr [0.95 – 1.05] p.u.
3
Buses: 18, 25 and 53
17
Branches: 19, 20, 31, 35, 36, 37, 41, 46, 54, 58, 59, 65, 66, 71, 73, 76 and 80 – 1250.8 MW, 336.4 MVAr [0.94 – 1.06] p.u.
33 – 50
14
9
130 – 64
Branches: 8, 32, 36, 51, 93, 95, 102, 107 and 127
– 4242 MW, 1439 MVAr [0.94 – 1.06] p.u.
3.1.2. Case 2: Minimization of cost considering multi-fuels Thermal generating plants have multi-fuel sources such as coal, natural gas and oil for different power output ranges of a unit. In such cases, the fuel cost function or curve is divided into piecewise quadratic cost functions (shown in Fig. 1) depending on the number and nature of fuels used. The cost considering multi-fuels for 𝑖th generator is mathematically expressed as: 𝑓𝑖 (𝑥, 𝑢) = 𝑎𝑖𝑘 + 𝑏𝑖𝑘 𝑃𝐺𝑖 + 𝑐𝑖𝑘 𝑃𝐺2 for 𝑓 𝑢𝑒𝑙𝑘 𝑖
(14)
within power output range 𝑃𝐺𝑚𝑖𝑛 ≤ 𝑃𝐺𝑖 ≤ 𝑃𝐺𝑚𝑎𝑥 ; 𝑘 being the fuel option. 𝑖𝑘 𝑖𝑘 The objective function reflecting the total cost is given by: ) (𝑁𝐺 ∑ (15) 𝑓𝑖 (𝑥, 𝑢) . 𝑓 (𝑥, 𝑢) = 𝑖=1
For the system under study, first two generating units are considered to have multi-fuel options. The corresponding coefficients and power output ranges (MW) are given in Table A.3 in Appendix. Costs for remaining 4-units are same as in case 1. 3.1.3. Case 3: Enhancement of voltage stability of the network Voltage stability problems are receiving growing attention in power systems as network collapses have been experienced in past due to voltage instability. Under normal condition and after being subjected to disturbance, the stability of a power system is characterized by its ability to maintain all bus voltages within acceptable limits. A system enters into a state of voltage instability when a disturbance, increase in load demand or change in system condition causes a progressive and uncontrollable decrease in voltage (Bhattacharya and Chattopadhyay, 2011). Systems with long transmission lines and heavy loading are more prone to voltage instability problem. In power system, improving voltage stability of a system is an important aspect. 𝐿-index of each bus serves as a good indicator of power system stability (Kessel and Glavitsch, 1986). The value of the index varies from 0 to 1, with 0 being the no load case while 1 signifies voltage collapse. If a power system has NL number of load (PQ) buses and NG number of generator (PV) buses, the value of L-index 𝐿𝑗 of bus 𝑗 is defined as: | | 𝑁𝐺 ∑ | 𝑉 | 𝐿𝑗 = ||1 − 𝐹𝑗𝑖 𝑖 || , where 𝑗 = 1, 2, … , 𝑁𝐿 𝑉𝑗 | | 𝑖=1 | | [ ]−1 and 𝐹𝑗𝑖 = − 𝑌𝐿𝐿 [𝑌𝐿𝐺 ]
Fig. 1. Comparison of cost curves of a generator with single fuel and with multi-fuels.
in (17). ][ ] [ ] [ 𝐼𝐿 𝑌 𝑌 𝑉𝐿 = 𝐿𝐿 𝐿𝐺 . 𝐼𝐺 𝑌𝐺𝐿 𝑌𝐺𝐿 𝑉𝐺
(17)
The L-index is calculated for all load buses and maximum value out of those acts as the global indicator for the system stability. Therefore, the objective function of system stability is given by: (18)
𝑓 (𝑥, 𝑢) = 𝐿𝑚𝑎𝑥 = max(𝐿𝑗 ), 𝑤ℎ𝑒𝑟𝑒 𝑗 = 1, 2, …, 𝑁𝐿.
3.1.4. Case 4: Minimization of emission Generation of electrical power from conventional sources of energy emits harmful gases into the environment. The amount of SOx , NOx emission in tonnes per hr (t/h) increases with increase in generated power (in p.u. MW) following the relationship given in Eq. (19). Minimization of emission is set as the objective of OPF.
(16)
where, sub-matrices 𝑌𝐿𝐿 and 𝑌𝐿𝐺 are obtained from system YBUS matrix after separating load (PQ) buses and generator (PV) buses as represented
𝑓 (𝑥, 𝑢) = Emission =
𝑁𝐺 [ ∑ 𝑖=1
84
(𝛼𝑖 + 𝛽𝑖 𝑃𝐺𝑖 + 𝛾𝑖 𝑃𝐺2 ) × 0.01 + 𝜔𝑖 𝑒(𝜇𝑖 𝑃𝐺𝑖 ) 𝑖
]
(19)
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
3.1.8. Case 8: Minimization of fuel cost and voltage deviation Voltage deviation is a measure of voltage quality in the network. The index of deviation is also important from security aspect. The indicator is formulated as cumulative deviation of voltages of all load buses (PQ buses) in the network from nominal value of unity. Mathematically it is expressed as: (𝑁𝐿 ) ∑| | 𝑉𝐷 = (23) |𝑉𝐿𝑝 − 1| . | | 𝑝=1
The combined objective function of fuel cost and voltage deviation is: 𝑓 (𝑥, 𝑢) =
(𝑁𝐺 ∑ 𝑖=1
) 𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2 𝑖
+ 𝜆𝑉 𝐷 × 𝑉 𝐷
(24)
where, weight factor 𝜆𝑉 𝐷 is assigned a value of 100 as in Bouchekara et al. (2016b) and Mohamed et al. (2017). 3.1.9. Case 9: Minimization of fuel cost and enhancement of voltage stability This objective function aims to minimize both fuel cost and enhance voltage stability of the system. The multiple objectives are converted to single objective as: ) (𝑁𝐺 ∑ 𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2 + 𝜆𝐿 × 𝐿𝑚𝑎𝑥 (25) 𝑓 (𝑥, 𝑢) =
Fig. 2. Cost curves of a generator with and without valve point effect.
𝑖=1
𝑖
where, 𝛼𝑖 , 𝛽𝑖 , 𝛾𝑖 , 𝜔𝑖 and 𝜇𝑖 are all emission coefficients provided in Table A.2 in Appendix for 30-bus system.
where, 𝐿𝑚𝑎𝑥 is calculated using Eq. (18) and selected value of weight factor 𝜆𝐿 is 100 (Mohamed et al., 2017).
3.1.5. Case 5: Minimization of real power loss The power loss in transmission system is unavoidable as the lines have inherent resistance. The real power loss (in MW) to be minimized is expressed as:
3.1.10. Case 10: Minimization of fuel cost, emission, voltage deviation and losses Four objectives are combined in this case study. Fuel cost, emission, voltage deviation and real power loss in the network are all minimized simultaneously. The objective function is given by: ) (𝑁𝐺 ∑ 𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2 𝑓 (𝑥, 𝑢) =
𝑓 (𝑥, 𝑢) = 𝑃𝑙𝑜𝑠𝑠 =
𝑛𝑙 ∑
[ ] 𝐺𝑞(𝑖𝑗) 𝑉𝑖2 + 𝑉𝑗2 − 2𝑉𝑖 𝑉𝑗 cos(𝛿𝑖𝑗 )
(20)
𝑞=1
where, 𝛿𝑖𝑗 = 𝛿𝑖 − 𝛿𝑗 , is the difference in voltage angles between bus 𝑖 and bus 𝑗 and 𝐺𝑞(𝑖𝑗) is the transfer conductance of branch 𝑞 connecting buses 𝑖 and 𝑗.
𝑖=1
( ( ))| | 𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2 + ||𝑑𝑖 × sin 𝑒𝑖 × 𝑃𝐺𝑚𝑖𝑛 − 𝑃𝐺𝑖 || 𝑖 𝑖 | | 𝑖=1
𝑁𝐺 ∑
3.2. IEEE 57-bus system General system data of 57-bus system are provided in Table 1. A total of four case studies are performed for the system with half of those being single objective and rest half are multi-objective optimizations. 3.2.1. Case 11: Minimization of fuel cost This basic objective function aims to minimize fuel cost for 57-bus system with optimal power flow. The equation of objective is same as (13). The values of generator cost (Zimmerman et al., 0000) and emission coefficients (Chaib et al., 2016), used in the calculation of cost and emission respectively, are provided in Table A.4 in Appendix.
(21)
where, 𝑑𝑖 and 𝑒𝑖 are the coefficients that represent the valve-point loading effect. The factors used for calculations are provided in Table A.2 in Appendix.
3.2.2. Case 12: Minimization of fuel cost and voltage deviation The goal of the objective function is to minimize both fuel cost and voltage deviation simultaneously. The converted single objective function follows Eq. (24) with weight factor 𝜆𝑉 𝐷 is assigned a value of 100, same as in the case for 30-bus system.
3.1.7. Case 7: Minimization of fuel cost and real power loss The multi-objective case of minimizing fuel cost and real power loss is converted to single objective by multiplying a weight factor to one of the objectives. It is represented as: 𝑓 (𝑥, 𝑢) =
𝑁𝐺 ∑ 𝑖=1
𝑎𝑖 + 𝑏𝑖 𝑃𝐺𝑖 + 𝑐𝑖 𝑃𝐺2 + 𝜆𝑃 × 𝑃𝑙𝑜𝑠𝑠 𝑖
(26)
The weight factors are selected as in Mohamed et al. (2017) with 𝜆𝐸 = 19, 𝜆𝑉 𝐷 = 21 and 𝜆𝑃 = 22 to balance among the objectives.
3.1.6. Case 6: Minimization of fuel cost considering valve point effect Valve-point effect needs to be considered for more realistic and precise modeling of fuel cost function. The generating units with multi-valve steam turbines exhibit a greater variation in the fuel-cost functions (Bouchekara et al., 2016b). The valve loading effect of multivalve steam turbines is modeled as sinusoidal function, the absolute value of which is added to the basic cost function. The real cost curve function of the steam plant becomes non-continuous as in Fig. 2. The objective of minimization of generating fuel cost with valvepoint effect is given by Biswas et al. (2017): 𝑓 (𝑥, 𝑢) =
𝑖
+ 𝜆𝐸 × 𝐸𝑚𝑖𝑠𝑠𝑖𝑜𝑛 + 𝜆𝑉 𝐷 × 𝑉 𝐷 + 𝜆𝑃 × 𝑃𝑙𝑜𝑠𝑠 .
3.2.3. Case 13: Minimization of fuel cost and enhancement of voltage stability The formulation of the objective function, consisting of both fuel cost and voltage stability, is same as in Eq. (25). The selected weight factor 𝜆𝐿 is also 100.
(22)
where, 𝑃𝑙𝑜𝑠𝑠 is the real power loss in the network calculated by using Eq. (20) and value of factor 𝜆𝑃 is chosen as 40, same as in Mohamed et al. (2017). 85
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
∙ 𝑋𝑖 is feasible but 𝑋𝑗 is infeasible ∙ 𝑋𝑖 and 𝑋𝑗 are both feasible and 𝑋𝑖 has better fitness value than 𝑋𝑗 ∙ 𝑋𝑖 and 𝑋𝑗 are both infeasible, but 𝑋𝑖 has a smaller overall constraint violation (computed by Eq. (29))
3.2.4. Case 14: Minimization of voltage deviation The objective function of minimizing cumulative voltage deviation of load buses from nominal value of 1 p.u. follows Eq. (23). 3.3. IEEE 118-bus system
Therefore, feasible individuals are always considered better than the infeasible individuals in this technique. Two feasible solutions are compared based on their objective function values only while two infeasible solutions are compared based on their overall constraint violations only. Comparison of infeasible solutions based on the overall constraint violation targets to push the infeasible solutions towards feasible region. Comparison of two feasible solutions on the objective value facilitates the convergence to optimal solution.
General system data of 118-bus system are provided in Table 1. A couple of single objective study cases are performed for the system. 3.3.1. Case 15: Minimization of fuel cost This basic objective function of minimizing fuel cost follows Eq. (13). The values of generator cost coefficients are derived from Zimmerman et al. (0000). 3.3.2. Case 16: Minimization of real power loss The objective aims to minimize real power loss in the 118-bus network. Power loss is calculated using Eq. (20). Network data are all taken from Zimmerman et al. (0000).
(b) Self-adaptive penalty (SP) Self-adaptive penalty (Farmani and Wright, 2003; Tessema and Yen, 2006) is an improved version of static penalty. In self-adaptive penalty, two types of penalties are added to each infeasible individual to identify the best infeasible individuals in the current population. The magnitudes of penalties added are controlled by the number of feasible individuals currently present in the combined population. If the number of feasible individuals is less, a higher amount of penalty is added to infeasible individuals with a higher amount of constraint violation. On contrary, if the number of feasible individuals is significantly large, then infeasible individuals with high fitness values will have small penalties added to their fitness values. These two penalties will allow the algorithm to switch between finding more feasible solutions and searching for the optimum solution at any time during the search process. Unlike static penalty, self-adaptive penalty requires no parameter tuning. The final fitness value based on which the population members are ranked is given as 𝐹 (𝑋) = 𝑑 (𝑋) + 𝑝(𝑋), where 𝑑 (𝑋) is the distance value and 𝑝 (𝑋) is the penalty value. The distance value is computed as follows: { 𝜈 (𝑋) , if 𝑟𝑓 = 0 (30) 𝑑 (𝑋) = √ 2 ′′ 2 [𝑓 (𝑋)] + [𝜈 (𝑋)] , otherwise 𝑓 (𝑋) − 𝑓𝑚𝑖𝑛 𝑓 ′′ (𝑋) = (31) 𝑓𝑚𝑎𝑥 − 𝑓𝑚𝑖𝑛 Number of feasible individuals (32) 𝑟𝑓 = Population size
4. Various constraint handling (CH) techniques integrated with DE As evolutionary algorithms (EAs) are stochastic in nature, infeasible solutions should not be straightaway discarded. A proper constraint handling method when used in conjunction with compatible evolutionary algorithm, it can drive the search process towards global feasible optima by making use of the information present in the infeasible solutions. Among the different constraint handling methods, the most commonly employed method is static penalty, where a penalty is added to the fitness of an infeasible solution for violating the constraints. Despite of its simplicity and ease of implementation, the performance of static penalty method heavily depends on the penalty factor which needs to be fine-tuned by trial and error. Therefore, few constraint handling methods that can be employed with EAs have been proposed and are summarized below. A constrained optimization problem with 𝑛 parameters to be optimized is usually written as a nonlinear programming problem of the following form (Qin et al., 2009; Mallipeddi and Suganthan, 2010): ( ) Minimize ∶ 𝑓 (𝑋) , 𝑋 = 𝑥1 , 𝑥2 , … , 𝑥𝑛 and 𝑋 ∈ 𝑆 subject to ∶
(27)
𝑔𝑖 (𝑋) ≤ 0, 𝑖 = 1, … , 𝑝 ℎ𝑗 (𝑋) = 0, 𝑗 = 𝑝 + 1, … , 𝑚
where, 𝜈 (𝑋) is the overall constrain violation as defined in Eq. (29), 𝑓𝑚𝑎𝑥 and 𝑓𝑚𝑖𝑛 are the maximum and minimum values of the objective function 𝑓 (𝑋) in the current combined population. The penalty value is ( ) defined as 𝑝 (𝑋) = 1 − 𝑟𝑓 𝑀 (𝑥) + 𝑟𝑓 𝑁(𝑥) where { 0, if 𝑟𝑓 = 0 𝑀 (𝑋) = and 𝜈 (𝑋) , otherwise { 0, if 𝑋 is a feasible individual 𝑁 (𝑋) = 𝑓 ′′ (𝑋) , if 𝑋is an infeasible individual.
where, 𝑆 is the entire search space, 𝑝 is the number of inequality constraints and (𝑚−𝑝) is the number of equality constraints. The equality constraints are transformed into inequality constraints and thus total constraint is represented as: [ ] { max g𝑖 (𝑋) , 0 , 𝑖 = 1, … , 𝑝 ] [ 𝐺𝑖 (𝑋) = (28) | | max |ℎ𝑖 (𝑋)| − 𝛿, 0 , 𝑖 = 𝑝 + 1, … , 𝑚 where, 𝛿 is a tolerance parameter for the equality constraints. The fitness function 𝑓 (𝑋) is to be minimized so that the optimal solution obtained in the process satisfies the inequality constraint 𝐺𝑖 (𝑋). The overall constraint violation for an infeasible solution is the weighted mean of all constraints and it is defined as: ∑𝑚 𝑤𝑖 (𝐺𝑖 (𝑋)) 𝜈(𝑋) = 𝑖=1∑𝑚 (29) 𝑖=1 𝑤𝑖
(c) Ensemble of constraint handling techniques (ECHT) Generally, each constrained optimization problem is unique in terms of the ratio between feasible search space and the whole search space, multi-modality and the nature of constraint functions. While solving a particular problem, the evolution paths can be different during different runs even if the same algorithm is used because of the stochastic nature of the EAs. In other words, the search process passes through different phases at different points during the evolution. Therefore, depending on several factors such as the ratio between feasible search space and the whole search space, multi-modality of the problem, nature of equality/inequality constraints, the chosen EA and global exploration/local exploitation stages of the search algorithm, different constraint handling methods can be effective during different stages of the search process. Due to the strong interactions among these diverse factors and the stochastic nature of the evolutionary algorithms, it is not straightforward to determine which constraint handling method is the best during a particular stage of the evolution to solve a given
where, 𝑤𝑖 (= 1∕𝐺𝑚𝑎𝑥,𝑖 ) is a weight parameter, 𝐺𝑚𝑎𝑥,𝑖 is the maximum value for violation of constraint 𝐺𝑖 (𝑋) recorded thus far. In the algorithm, 𝑤𝑖 is set as 1∕𝐺𝑚𝑎𝑥,𝑖 , which varies during the evolution to balance the contribution of each constraint in the problem regardless of the different numerical ranges of all constraints. A brief overview of the two constraint handling techniques used for OPF problem is presented herein. (a) Superiority of feasible solutions (SF) In SF (Deb, 2000), solution 𝑋𝑖 is considered superior to solution 𝑋𝑗 when: 86
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
the other constraint handling methods. Therefore, an offspring produced by a particular constraint handling method may be rejected by its own population, but could be accepted by the populations of other constraint handling methods. Hence, in ECHT every function call is utilized effectively making it useful for objective/constraint functions that are computationally expensive. In ECHT, if a particular constraint handling technique is best suited for a given problem and the search method at a particular point in the search process, the offspring population produced by the population of that constraint handling method will dominate the other populations and enter other populations to become parents in the subsequent generations. Therefore, ECHT takes the onus of choosing the best constraint handling technique and tuning the associated parameter values for a particular problem. A summary of steps involved in ECHT-DE is provided in Table 2. Flowchart with basic steps of the algorithm is also presented in Fig. 3. Dimension 𝐷 of OPF problem is 24, 33 and 130 for IEEE 30-bus, 57-bus and 118-bus systems, respectively. Population size 𝑁𝑃 and maximum number of function evaluation, 𝑚𝑎𝑥𝑒𝑣𝑎𝑙 are varied for different cases and those details are mentioned in Section 5. The proposed algorithms are developed using MATLAB software and simulations are carried out on a computer with Intel Core i5 CPU @2.7 GHz and 4GB RAM. Simulation results are discussed in subsequent section.
problem using a given EA. Motivated by these observations, ECHT was developed (Mallipeddi and Suganthan, 2010) to implicitly benefit from the match between constraint handling methods, characteristics of the problem being solved, chosen EA and the exploration–exploitation stages of the search process. (d) Differential evolution (DE) Differential Evolution (DE), introduced by Storn and Price (Storn and Price, 1997), is a stochastic, population based optimization algorithm. The algorithm employs crossover and mutation operators that help the individuals in the population to evolve and improve their fitness. These individuals are evaluated and those that perform better are selected to compose the population in the next generation. A brief description of the steps involved in DE algorithm is presented in this section. Initialization DE optimization process initializes the population of candidate solutions with assignment of random values to each decision vector in the population. These values must be within the feasible bounds (between maximum & minimum) of the decision variables. Initialization of 𝑗th component of 𝑘th decision vector is done as: 𝑋𝑘,0 (𝑗) = 𝑥𝑚𝑖𝑛 (𝑗) + 𝑟𝑎𝑛𝑑𝑘,𝑗 [0, 1] × [𝑥𝑚𝑎𝑥 (𝑗) − 𝑥𝑚𝑖𝑛 (𝑗)]
(33) 5. Results and comparison
where, 𝑟𝑎𝑛𝑑𝑘,𝑗 [0, 1] is a random number lying between 0 and 1; 𝑗 = 1, 2, … , 𝐷 where 𝐷 is the dimension of the decision vector.
Simulation results of application of SF-DE, SP-DE and ECHT-DE algorithms are tabulated, and also detail system by system explanation and comparative analysis are provided in this section. Table 3 summarizes the objectives of various case studies performed for the test systems with all the 3 algorithms. The parameters to be optimized for a specific case are ticked in corresponding cells of the table. A population size (𝑁𝑃 ) of 30 is considered for all the study cases of the systems under study. The population size is selected after several trial runs. Increasing number of population members adversely affects the algorithm in finding feasible solutions. Again, for single objective and multi-objective optimizations of IEEE 30-bus system, 30 000 function evaluations (𝑚𝑎𝑥𝑒𝑣𝑎𝑙 for system 1) are performed in a single complete run of the algorithm. A total of 42 000 function evaluations (𝑚𝑎𝑥𝑒𝑣𝑎𝑙 for system 2) of objective functions are carried out for IEEE 57-bus system in one run. Due to large number of control variables in 118-bus system, 300 000 function evaluations (𝑚𝑎𝑥𝑒𝑣𝑎𝑙 for system 3) are set as the stopping criterion in a study case. Each case of 30-bus and 57-bus systems (i.e. case 1 to case 14) is run 30 times independently to perform statistical analysis and evaluate the performance of the 3-constraint handling techniques.
Mutation In second step, mutation operator of DE creates a donor/mutant vector 𝑉𝑘,𝐺 corresponding to each population member or target vector 𝑋𝑘,𝐺 in the current generation (the subscript ‘G’ denotes parameter at G-th generation). The strategy for mutation used here is ‘DE/rand/1’: 𝑉𝑘,𝐺 = 𝑋𝑅𝑘1,𝐺 + 𝐹 (𝑋𝑅𝑘2,𝐺 − 𝑋𝑅𝑘3,𝐺 ).
(34)
The indices 𝑅𝑘1, 𝑅𝑘2 and 𝑅𝑘3 are mutually exclusive and chosen randomly from the population range. The scaling (mutation) factor 𝐹 is a positive control parameter for scaling the difference vectors. Crossover The donor vector mixes its components with the target vector 𝑋𝑘,𝐺 to form the trial/offspring vector 𝑈𝑘,𝐺 = [𝑈𝑘,𝐺 (1) , 𝑈𝑘,𝐺 (2) , … .., 𝑈𝑘,𝐺 (𝐷)] through crossover. Binomial crossover is adopted here and the scheme operates on each variable whenever a randomly generated number between 0 and 1 is less than or equal to a pre-fixed crossover rate 𝐶𝑅. For an element, the scheme is expressed as: { 𝑉𝑘,𝐺 (𝑗) if 𝑗 = 𝑗𝑟𝑎𝑛𝑑 or 𝑟𝑎𝑛𝑑𝑘,𝑗 [0, 1] ≤ 𝐶𝑅, 𝑈𝑘,𝐺 (𝑗) = (35) 𝑋𝑘,𝐺 (𝑗) otherwise, where 𝑗 = 1, … , 𝐷
5.1. Statistical comparison among ECHT-DE, SF-DE and SP-DE
where, 𝑗𝑟𝑎𝑛𝑑 is a randomly chosen natural number in {1, 2, . . . , D }. The new offspring vector, 𝑈𝑘,𝐺 is evaluated and compared against parent vector 𝑋𝑘,𝐺 in terms of its fitness value and constraint violation. In case of single constraint handling method applied to a problem, the ordering of feasible solutions and operation on infeasible solutions among offspring and parent vectors are dictated by the rules set up by that constraint handling method, i.e. either by SF or SP as described in clause 4(a) and (b). In real-world problem such as OPF where a single run to find optimum settings of control variables with metaheuristics can take considerable time and resource, finding a suitable constraint handling method by trial-and-error may become difficult. The computation time needed to find an appropriate constraint handling method can be saved by using the ECHT method where two constraint handling techniques discussed in Section 4(a) and (b) and DE (Storn and Price, 1997) as the basic search algorithm are used. In ECHT, each constraint handling technique has its own population and parameters that produce offsprings for evaluation. The parent population corresponding to a particular constraint handling method not only competes with its own offspring population but also with offspring populations of
The statistical summary of 30 runs for each study case performed (i.e. case 1 to case 14) using different constraint handling techniques with DE as the basic search algorithm is presented in Table 4. The columns indicate the best, worst, mean and standard deviation values for the objective function in each case. It is evident that no single method is able to provide best fitness or best mean values in all the cases. Further, we carry out Wilcoxon signed rank test to compare performance between any two algorithms. The procedure applied for comparative performance analysis between algorithm ‘A’ and algorithm ‘B’ on a specific study case using signed rank test (Derrac et al., 2011) is described below. ➢ Gather all fitness values over 30 runs of the objective in a case study for both the algorithms ➢ Compute R+, the sum of ranks for runs in which algorithm ‘A’ outperforms algorithm ‘B’ ➢ Compute R−, the sum of ranks for runs in which algorithm ‘B’ outperforms algorithm ‘A’ 87
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 2 Steps in ECHT-DE. ECHT-DE INPUT
∙ Dimension of the problem, 𝐷 ∙ Number of population for each constraint handling method, 𝑁𝑃 ∙ Stopping criteria, maximum number of function evaluation, 𝑚𝑎𝑥𝑒𝑣𝑎𝑙 ∙ Maximum and minimum values of 𝐷-decision variables, in vector form 𝑋𝑚𝑎𝑥 and 𝑋𝑚in 1 D 1 D ∙ 𝑋𝑚𝑎𝑥 = [𝑋𝑚𝑎𝑥 , … , 𝑋𝑚𝑎𝑥 ] and 𝑋𝑚in = [𝑋min , … , 𝑋min ]
INITIALIZATION
∙ Set generation counter, 𝐺 = 0 and function evaluation counter 𝑛𝑓 𝑒𝑣𝑎𝑙 = 0. ∙ Set DE parameters: Mutation factor, 𝐹 = 0.7; Crossover rate, 𝐶𝑅 = 0.9 ∙ POP1: Create population 1 of 𝑁𝑃 individuals uniformly distributed between [𝑋𝑚𝑎𝑥 , 𝑋𝑚𝑖𝑛 ] ∙ POP2: Create population 2 of 𝑁𝑃 individuals uniformly distributed between [𝑋𝑚𝑎𝑥 , 𝑋𝑚𝑖𝑛 ] ∙ Each of the two constraint handling techniques (SF, SP) has its own population POP1 and POP2 ∙ Initialize parameters PAR1 and PAR2 for SF and SP respectively according to Section 4(a) and (b)
EVALUATION
∙ Evaluate objective function, constraint function and constraint violation using Eq. (27) to Eq. (29) for each individual 𝑋1k , ∀𝑘 ∈ {1, … , 𝑁𝑃 } of POP1 ∙ Evaluate objective function, constraint function and constraint violation using Eq. (27) to Eq. (29) for each individual 𝑋2k , ∀𝑘 ∈ {1, … , 𝑁𝑃 } of POP2 ∙ Increase function evaluation counter 𝑛𝑓 𝑒𝑣𝑎𝑙 by 2 × 𝑁𝑃
ALGORITHM LOOP: STEP 1
∙ Update parameter values PAR1 for SF constraint handling methods following Section 4(a) ∙ Update parameter values PAR2 for SP constraint handling methods following Section 4(b)
STEP 2
∙ Perform mutation Eq. (34) and crossover Eq. (35) for each individual of POP1; generate offspring OFS1 ∙ Perform mutation Eq. (34) and crossover Eq. (35) for each individual of POP2; generate offspring OFS2 ∙ Sample mutation and crossover strategy for POP1 at generation 𝐺: Mutant vector, 𝑉 1𝑘,𝐺 = 𝑋1𝑅𝑘1,𝐺 + 𝐹 (𝑋1𝑅𝑘2,𝐺 − 𝑋1𝑅𝑘3,𝐺 ), 𝑘 ∈ {1, … , 𝑁𝑃 } { 𝑉 1𝑘,𝐺 (𝑗) if 𝑗 = 𝑗𝑟𝑎𝑛𝑑 or 𝑟𝑎𝑛𝑑𝑘,𝑗 [0, 1] ≤ 𝐶𝑅, Trial vector element, 𝑈 1𝑘,𝐺 (𝑗) = 𝑋1𝑘,𝐺 (𝑗) otherwise, where 𝑗 = 1, … , 𝐷 where 𝑅𝑘1, 𝑅𝑘2 and 𝑅𝑘3 are mutually exclusive integers randomly chosen from population range {1, … , 𝑁𝑃 }; 𝑗𝑟𝑎𝑛𝑑 is a randomly chosen natural number within {1, … , 𝐷}
STEP 3
∙ Compute objective function, constraint function and constraint violation using Eq. (27) to Eq. (29) for each individual 𝑋1′k , ∀𝑘 ∈ {1, … , 𝑁𝑃 } of OFS1 ∙ Compute objective function, constraint function and constraint violation using Eq. (27) to Eq. (29) for each individual 𝑋1′k , ∀𝑘 ∈ {1, … , 𝑁𝑃 } of OFS2 ∙ Increase function evaluation counter 𝑛𝑓 𝑒𝑣𝑎𝑙 by 2 × 𝑁𝑃
STEP 4
∙ Combine parent population with all offsprings to create groups. ∙ Group1: (POP1+OFS1+OFS2) and Group2: (POP2+OFS1+OFS2)
STEP 5
∙ In selection step, parent populations POP1 and POP2 for the next generation are selected from Group1 and Group2 respectively. In a Group (say Group1), since OFS1 is produced by POP1 by mutation and crossover, DE’s selection based on competition between parent and its offspring is employed when POP1 competes with OFS1. But when POP1 competes with OFS2, produced by other populations, each member in POP1 competes with a randomly selected offspring from OFS2.
STEP 6
∙ Is stopping criterion 𝑚𝑎𝑥𝑒𝑣𝑎𝑙 reached? ∙ If yes, STOP ∙ If not, increase generation counter by 1, i.e. 𝐺 = 𝐺 + 1. Go to algorithm loop STEP 1.
Table 3 Summary of case studies. Case no.
IEEE 30-bus system Basic fuel cost
Case 1 Case 2 Case 3 Case 4 Case 5 Case 6 Case 7 Case 8 Case 9 Case 10 Case 11 Case 12 Case 13 Case 14 Case 15 Case 16
Multi-fuel cost
IEEE 57-bus system Voltage stability
Emission
Power loss
Voltage deviation
Valve-point effect
Basic fuel cost
Voltage stability
✔ ✔ ✔
✔
IEEE 118-bus system Voltage deviation
Basic fuel cost
Power loss
✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔
✔ ✔ ✔ ✔ ✔
✔
✔ ✔ ✔ ✔ ✔
➢ Calculate 𝑝-value which determines the significance of results in a statistical hypothesis test. The smaller the 𝑝-value, the stronger is the proof against null hypothesis 𝐻0.
8. In the remaining cases, SF-DE is found to be superior. ECHT-DE outperforms SP-DE in cases 1, 3, 6, 7 and 9 for 30-bus system. However, for all cases of 57-bus system, SP-DE outdoes ECHT-DE. An interesting point to note is in the performance comparison between SF-DE and SP-DE. While in cases 1, 3, 6, 7 and 9, SP-DE is inferior to SF-DE, in all remaining cases it is equal or better. This very fact reiterates the need and justification for development of ensemble method. The factors influencing the performance of a CH technique are the nature of the objective function, constraints, the feasible and infeasible regions for the
The results obtained using Wilcoxon signed rank test are presented in Table 5. The column 𝐻0 defines whether null hypothesis is valid or not. If null hypothesis is valid (i.e. 𝐻0 = ‘Yes’ with a significance level, 𝛼 = 0.05), the performance of two methods is statistically same for the study case. ECHT-DE and SF-DE perform equally in cases 4 and 88
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Fig. 3. Flowchart for implementation of ECHT-DE.
Table 4 Statistical summary of case studies using 3-CH techniques. Case no.
ECHT-DE
SF-DE
SP-DE
Best
Worst
Mean
Std dev
Best
Worst
Mean
Std dev
Best
Worst
Mean
Case 1
800.4148
800.4258
800.4206
0.0026
800.4131
800.4192
800.4151
0.0015
800.4293
800.4684
800.4413
0.0100
Case 2
646.4532
646.5423
646.4958
0.0232
646.4390
646.5117
646.4776
0.0178
646.4314
646.496
646.4565
0.0169
Case 3
0.13632
0.13721
0.13693
0.0002
0.13671
0.13704
0.13685
0.0001
0.13743
0.13857
0.13777
0.0002
Case 4
0.20482
0.20482
0.20482
0
0.20482
0.20482
0.20482
0
0.20482
0.20482
0.20482
0
Case 5
3.0850
3.0871
3.0858
0.0005
3.0845
3.0857
3.0849
0.0003
3.0844
3.08536
3.0848
0.0003
Case 6
832.1356
832.2239
832.1811
0.0222
832.0882
832.1291
832.1056
0.0105
832.4813
832.876
832.655
0.0963
Case 7
1040.151
1040.233
1040.181
0.0213
1040.125
1040.162
1040.14
0.0096
1040.134
1040.337
1040.239
0.0444
Case 8
813.1742
813.4095
813.2470
0.049
813.1956
813.3376
813.2585
0.0361
813.1959
813.2643
813.2306
0.0181
Case 9
814.1708
814.2001
814.1843
0.0075
814.1649
814.1956
814.1767
0.0063
814.1841
814.2273
814.2017
0.0121
Case 10
964.1331
964.1564
964.1437
0.0061
964.1254
964.1418
964.1307
0.0038
964.1234
964.1399
964.1276
0.0034
Case 11
41670.56
41680.09
41675.33
2.2866
41667.85
41684.29
41671.16
3.5896
41667.82
41678.96
41670.44
2.2474
Case 12
41776.48
41783.22
41778.87
1.7681
41775.09
41778.39
41776.16
0.6806
41774.75
41776.19
41775.3
0.3724
Case 13
41699.25
41727.03
41706.5
6.1009
41695.55
41712.35
41699.51
3.7249
41696.54
41700.97
41698.17
1.0153
Case 14
0.60416
0.64184
0.62012
0.0094
0.59584
0.61677
0.60602
0.0059
0.59267
0.60284
0.59800
0.0028
solutions, and the way the CH technique works. In SP, higher amount of penalty is added to infeasible individuals if a few feasible individuals are present in the population (Tessema and Yen, 2006). On the other hand, basic SF technique has the problem of maintaining population diversity and may stagnate especially when ratio between feasible region and search space is small and initial population has a few feasible solutions (Coello, 2002). The limited number of feasible solutions during initial phases of search process can deter SF from attaining
Std dev
optimum values in the trial runs. In summary, it is a tedious task to select appropriate CH method for a practical problem evaluated under different scenarios. Gauging extent of non-linearity in the objective function or assessing nature of constraints or delimiting feasible region is not straightforward. Ensemble method has the advantage of accessing and operating on the population members of all CH techniques it is composed of. Due to probabilistic nature of DE, ECHT-DE is able to attain close to, if not the best fitness value in most trial runs. 89
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 5 Wilcoxon signed rank test results. Case no. Case 1 Case 2 Case 3 Case 4 Case 5 Case 6 Case 7 Case 8 Case 9 Case 10 Case 11 Case 12 Case 13 Case 14
ECHT-DE vs SF-DE
ECHT-DE vs SP-DE
SF-DE vs SP-DE
R+
R−
𝑝-value
𝐻0
R+
R−
𝑝-value
𝐻0
R+
R−
𝑝-value
H0
1 94 103.5 0 8.5 0 0 303 31 5 53 0 11 13.5
464 371 361.5 0 456.5 465 465 162 434 460 412 465 454 451.5
1.92× 10−6 4.39 × 10−3 0.00793 1 4.07 × 10−6 1.73 × 10−6 1.73 × 10−6 0.14704 3.41 × 10−5 2.88 × 10−6 2.22 × 10−4 1.73 × 10−6 5.22 × 10−6 6.65 × 10−6
No No No Yes No No No Yes No No No No No No
465 16 465 0 1 465 446 176 425 0 28 0 2 0
0 449 0 0 464 0 19 289 10 465 437 465 463 465
1.73 × 10−6 8.47 × 10−6 1.73 × 10−6 1 1.92 × 10−6 1.73 × 10−6 1.13 × 10−5 0.24519 7.23 × 10−6 1.73 × 10−6 2.60 × 10−5 1.73 × 10−6 2.13 × 10−6 1.73 × 10−6
No No No Yes No No No Yes No No No No No No
465 67 465 0 162.5 465 464 83 465 72 191 14 168 12
0 398 0 0 302.5 0 1 382 0 393 274 451 297 453
1.73 × 10−6 6.64 × 10−4 1.73 × 10−6 1 0.14986 1.73 × 10−6 1.92 × 10−6 2.10 × 10−3 1.73 × 10−6 9.63 × 10−4 0.39333 6.98 × 10−6 0.18462 5.75 × 10−6
No No No Yes Yes No No No No No Yes No Yes No
5.2. System 1: IEEE 30-bus test system In this sub-section initially in Table 6, we list down settings of all control and state (dependent) variables alongwith their allowable ranges for the best fitness value obtained for the objective function in a case study pertaining to 30-bus test system using one of the 3CH techniques (ECHT, SF and SP). Active power of swing generator (𝑃G1 for 30 and 57-bus systems) and reactive power of all generators are state or dependent variables, treated as constraints in the optimization. Values of these variables are listed here to show that the CH techniques duly comply with the limits of these constraining variables. Limits of allowable reactive power for generators are derived from Zimmerman et al. (0000). Last row of the table indicates the approximate time needed by the processor (CPU) for one run of a study case. Table 7 presents comparison of results obtained using the 3-CH techniques in present study with past studies for single objective functions. Table 8 provides similar comparison for multi-objective cases. For case 1 of optimizing basic fuel cost, ECHT-DE and SF-DE algorithms lead to fuel costs of 800.4148 $/h and 800.4131 $/h, respectively satisfying all system constraints, more precisely conforming to the important inequality constraints on generator reactive power, line capacity and load bus voltage. Among all inequality constraints, constraint on load bus voltage is found to be critical as the operating voltages of load buses are often found to be near the limits. Some of the recent literatures recorded better results than what have been achieved by the 3-methods in present study. A careful observation of those results reveals violation of voltage limits [0.95 p.u.−1.05 p.u.], thus rendering the solutions infeasible. Cursory glance at voltage profiles of load buses leads to the presumption that the violation has all likely occurred exceeding the upper limit. This situation of overvoltage is undesirable as it may stress the system and often lead to failure of the connected equipment. For the 30-bus system, a total of 24 load buses (PQ bus) are present. Theoretically, if all these buses operate at the limits, the maximum possible cumulative absolute value of voltage deviation (VD in Eq. (23)) would be 1.2 p.u. (i.e. 24×0.05 p.u.). The reported VD values in some literatures are higher than 1.2 p.u. in many cases. As it turns out, the contentious results are output of static penalty method which is sensitive to selection of penalty coefficient. In comparison Tables 7 and 8, such cases are highlighted with footnotes. In this specific problem of OPF, limits on all state or dependent variables are duly satisfied in our study. In static penalty function approach, the constraints on these variables are violated in many cases, possibly without cognizance of the programmer. Moreover, the superiority of a proper CH technique lies in its ability to find the best result by operating near to the limits but not violating those. Figs. 4 and 5 indicate the load bus voltages for various study cases performed for IEEE 30-bus system. The voltage profile for a study case corresponds to the control variables obtained using the best performing method on that case as listed in Table 6. Clearly voltage limits are all conformed for the best solution.
Fig. 4. IEEE 30-bus system – voltage profiles of load buses for best solutions of case 1 to case 5.
Fig. 5. IEEE 30-bus system – voltage profiles of load buses for best solutions of case 6 to case 10.
In case 2 of minimizing cost with multiple fuels, the present study with SP-DE finds a cost of 646.4314 $/h and turns out to be among the best values when compared with results achieved by other algorithms. In case 3 of voltage stability indicator, i.e. minimization of maximum ‘L-index’ (𝐿𝑚𝑎𝑥 ) of the system load buses, ECHT-DE gives best output of 0.13632, marginally better than the other equivalent algorithms. Objective of emission in tonnes per hour (t/h) and other useful calculated parameters obtained by ECHT-DE, SF-DE and SP-DE in case 4 are almost same with those reported by algorithms MSA (Mohamed 90
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 6 Simulation results of best solutions for system 1: IEEE 30-bus system. Parameters Algorithm 𝑃G2 (MW) 𝑃G5 (MW) 𝑃G8 (MW) 𝑃G11 (MW) 𝑃G13 (MW) 𝑉1 (p.u.) 𝑉2 (p.u.) 𝑉5 (p.u.) 𝑉8 (p.u.) 𝑉11 (p.u.) 𝑉13 (p.u.) 𝑄c10 (MVAr) 𝑄c12 (MVAr) 𝑄c15 (MVAr) 𝑄c17 (MVAr) 𝑄c20 (MVAr) 𝑄c21 (MVAr) 𝑄c23 (MVAr) 𝑄c24 (MVAr) 𝑄c29 (MVAr) 𝑇11 (p.u.) 𝑇12 (p.u.) 𝑇15 (p.u.) 𝑇36 (p.u.) Fuel cost ($/h) Emission (t/h) 𝑃loss (MW) VD (p.u.) L-index (max) 𝑃G1 (MW) 𝑄G1 (MVAr) 𝑄G2 (MVAr) 𝑄G5 (MVAr) 𝑄G8 (MVAr) 𝑄G11 (MVAr) 𝑄G13 (MVAr) CPU time (s)
Min
Max
20 15 10 10 12 0.95 0.95 0.95 0.95 0.95 0.95 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.90 0.90 0.90 0.90
80 50 35 30 40 1.10 1.10 1.10 1.10 1.10 1.10 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 1.10 1.10 1.10 1.10
50 −20 −20 −15 −15 −10 −15
200 150 60 62.5 48.7 40 44.7
Case 1
Case 2
Case 3
Case 4
Case 5
Case 6
Case 7
Case 8
Case 9
Case 10
SF-DE 48.7282 21.3535 21.1986 11.9028 12.0021 1.0832 1.064 1.0333 1.038 1.0808 1.0418 3.6349 2.6015 3.8678 4.9883 4.184 4.9792 3.1586 4.9816 2.3958 1.051 0.9254 0.9649 0.9727 800.4131 0.36652 9.0104 0.92097 0.13786 177.2252 2.9154 19.2342 26.0422 28.0047 21.0932 −6.0223 133.1
SP-DE 54.9987 23.9573 34.9706 18.6602 17.5408 1.0753 1.0612 1.0314 1.0406 1.077 1.051 4.3281 0.2885 3.812 4.1986 4.0378 4.7194 4.0527 4.7473 2.3624 1.0231 0.9575 0.9765 0.9767 646.4314 0.28351 6.7276 0.91253 0.13832 140 −2.4871 17.6727 23.9799 27.182 16.2243 0.9936 122.2
ECHT-DE 79.7793 49.718 34.1924 28.8158 13.733 1.0542 1.051 1.0682 1.0568 1.0997 1.0741 2.9321 0.3526 1.5598 1.2874 0.0082 0.0327 0.0021 0.2923 0.0232 1.0445 0.9001 1.0035 0.9646 917.5916 0.2252 4.5224 0.9110 0.13632 81.684 −19.9831 −19.9361 54.341 48.4862 27.0418 18.6036 136.5
ECHT-DE 67.5933 50.0 35.0 30.0 40.0 1.0622 1.0565 1.0373 1.0439 1.0691 1.0507 1.5203 3.8625 4.2725 4.9037 3.8882 4.9426 3.1486 4.9819 2.5522 1.0677 0.901 0.9905 0.9771 944.3782 0.20482 3.2179 0.89240 0.13844 64.0253 −5.6506 8.0392 21.9063 28.345 19.4601 1.5932 138.2
SP-DE 80.0 50.0 35.0 30.0 40.0 1.0613 1.0572 1.0378 1.0442 1.0685 1.055 3.3113 0.8576 4.1219 4.9901 4.0549 5.0 2.8962 4.9919 2.2836 1.0656 0.9071 0.9908 0.9763 967.5962 0.20726 3.0844 0.90359 0.13832 51.4968 −5.5594 7.2517 21.7516 27.7104 18.6164 4.7497 136.4
SF-DE 44.8982 18.5075 10.023 10.0013 12.0006 1.0846 1.062 1.0285 1.0348 1.084 1.0516 2.9496 3.7232 3.9523 4.8533 3.5919 4.9914 3.6448 5.0 2.2825 1.0068 0.9911 0.9796 0.9773 832.0882 0.43730 10.6387 0.84935 0.13934 198.6081 4.9696 16.9485 24.3502 28.5989 18.1039 1.3342 137.6
SF-DE 55.4313 38.1719 35.0 30.0 26.7077 1.0687 1.0581 1.0348 1.0427 1.0753 1.0507 2.8395 0.8649 4.3812 4.9981 3.9664 4.9924 2.7251 4.9889 2.2225 1.0634 0.9076 0.9811 0.9737 859.1458 0.2289 4.5245 0.92731 0.13785 102.6169 −3.3869 10.9563 23.0135 27.5277 20.7233 1.0117 137.9
ECHT-DE 48.6293 21.6905 22.6953 11.8749 12.0233 1.0401 1.024 1.0144 1.0059 1.0424 0.9881 4.9484 0.5117 4.9733 0.0085 4.9924 4.9014 4.9636 4.9872 2.5552 1.0624 0.9002 0.9393 0.97 803.7198 0.36384 9.8414 0.09454 0.14888 176.3282 −4.3297 14.3453 46.9979 40.3091 21.4086 −14.6404 123.3
SF-DE 48.786 21.212 21.5819 11.8277 12.0002 1.0831 1.0641 1.0329 1.0379 1.0928 1.0455 0.7465 0.301 4.517 4.7206 3.7417 4.9716 2.7587 4.9129 2.0903 1.0265 0.9426 0.9667 0.9697 800.4203 0.36592 8.9985 0.93846 0.13745 176.9908 2.5344 20.1818 25.7862 28.8048 22.6468 −3.2481 130.4
SP-DE 52.528 31.4679 35.0 26.6793 21.0893 1.0723 1.0589 1.0317 1.0403 1.0305 1.0124 3.6781 3.9279 4.0292 4.9399 4.9874 5.0 4.1832 5.0 2.5508 1.0882 0.9529 1.021 1.0043 830.2123 0.25294 5.5857 0.29615 0.14756 122.2219 −1.0989 13.1559 22.963 26.8471 15.8058 0.7012 128.8
et al., 2017) and ARCBBO (Kumar and Premalatha, 2015). Case 5 by SF-DE and SP-DE result in real power losses of 3.0845 MW and 3.0844 MW, respectively, comparable with the results summarized in the table. Valve-point effect is considered in case 6 to arrive at an increase in cost than in case 1 with final value being 832.0882 $/h, obtained by SFDE. In summary, although variation in performance is observed among 3-CH techniques, output results of one or more methods employed in our study are among the best reported results on the problem of OPF. Moreover, the objective of this paper is not to prove superiority on mere numerical results, but to demonstrate strict compliance with the system constraints by the CH techniques. Case 7 to case 10 are for multi-objective cases of OPF for 30bus system. In these case studies, fitness of the combined objective function is the major deciding factor in ranking the results out of various optimization algorithms. For a meaningful comparison, fitness value of other methods is calculated and produced here using the same weight factor of the respective objective function. In multi-objective cases, a change in weight factor e.g. higher weight factor on fuel cost in case 7 will bring down the cost but will increase power loss. Quick review of Table 8 reveals that either one or more algorithms among ECHT-DE, SF-DE and SP-DE achieve best fitness values in all the cases. Though best fitness is reported by SF-DE in case 7, an intermediate value fuel cost, one of the constituting objectives, is achieved. The other objective of real power loss is the least when assessed with other algorithms. In case 8 of minimizing cost and voltage deviation (𝑉 𝐷), 𝑉 𝐷 achieved by ECHT-DE is the lowest among all other comparable methods. Optimum values of both fitness and cost are the lowest in run of case 9 by SF-DE algorithm. Four objectives of cost, emission, real power loss and voltage deviation are simultaneously minimized in case 10. Alongwith fitness value, SP-DE arrives at lowest values of cost and loss when compared with MSA (Mohamed et al., 2017) and MFO (Mohamed et al., 2017).
Fig. 6. Comparative convergence of case 1 for IEEE 30-bus system.
Graphical comparisons of convergences of 3-CH techniques for case 1, case 2 and case 6 of fuel cost related objective functions are illustrated in Figs. 6, 7 and 11 respectively. The convergence speeds are not remarkably different among the algorithms. However, rapid and sudden convergence is observed for both ECHT-DE and SP-DE during initial phase of the search process. SF-DE converges to the optimal solution more steadily. For cases 3, 4 and 5, the convergences are indicated in Figs. 8–10 respectively. The convergence is a little more uneven in case 3 due to the nature of the objective function. In general, the algorithms converge to optimal solutions somewhere around 15 000 function evaluations. Convergences of two-objective cases are provided in Fig. 12–14. For clarity, convergence of only one algorithm attaining best fitness 91
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 7 Comparison of 3-CH techniques with previous studies for single objective cases of IEEE 30-bus system. Case #
Algorithm
Fuel cost ($/h)
Emission (t/h)
𝑃loss (MW)
VD (p.u.)
L-index (max)
Case 1
ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) AGSO (Daryani et al., 2016) ARCBBO (Kumar and Premalatha, 2015) SKH (Pulluri et al., 2016) BSA (Chaib et al., 2016) ICBO (Bouchekara et al., 2016b) APFPA (Mahdad and Srairi, 2016) FHSA (Pandiarajan and Babulal, 2016) GEM (Bouchekara et al., 2016a) DE (Shaheen et al., 2016) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) LTLBO (Ghasemi et al., 2015) GABC (Roy and Jadhav, 2015) BSA (Chaib et al., 2016) ICBO (Bouchekara et al., 2016b) ECHT-DE SF-DE SP-DE SKH (Pulluri et al., 2016) GEM (Bouchekara et al., 2016a) DE (Shaheen et al., 2016) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) DSA (Abaci and Yamacli, 2016) AGSO (Daryani et al., 2016) ARCBBO (Kumar and Premalatha, 2015) GEM (Bouchekara et al., 2016a) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) DSA (Abaci and Yamacli, 2016) ARCBBO (Kumar and Premalatha, 2015) APFPA (Mahdad and Srairi, 2016) GEM (Bouchekara et al., 2016a) ECHT-DE SF-DE SP-DE PSO (Bouchekara et al., 2016b) ICBO (Bouchekara et al., 2016b) BSA (Chaib et al., 2016) APFPA (Mahdad and Srairi, 2016)
800.4148 800.4131 800.4293 800.5099 801.75 800.5159 800.5141 799.0760a 799.0353a 798.9144a 799.914a 799.0463a 799.0827a 646.4532 646.4390 646.4314 646.8364 647.4315 647.03 646.1504a 645.1668a 917.5916 875.8929 842.015 814.0100 816.9095a 915.2172a 944.3782 944.3242 944.3323 944.5003 944.4086 953.629 945.1597 943.6358a 967.6001 967.6198 967.5962 967.6636 967.6493 967.6605 965.6590a 966.7473a 832.1356 832.0882 832.4813 832.6871 830.4531a 830.7779a 830.4065a
0.36593 0.36652 0.36806 0.36645 0.3703 0.3663 0.3662 0.3671 – – – 0.3665 – 0.28351 0.2835 0.28351 0.28352 0.2835 – 0.2833 – 0.2252 0.22801 0.24635 0.3740 0.2802 – 0.20482 0.20482 0.20482 0.20482 0.20583 0.2059 0.2048 0.2048a 0.20726 0.20726 0.20726 0.20727 0.20826 0.2073 – 0.2072 0.43765 0.43730 0.43651 – – 0.4377 –
8.9999 9.0104 9.0583 9.0345 – 9.0255 9.0282 8.6543 8.6132 8.5800 – 8.6257 8.63 6.722 6.7258 6.7276 6.8001 6.9347 6.8160 6.6233 6.3828 4.5224 4.6412 5.92 9.9056 6.2313 3.626 3.2179 3.2179 3.2178 3.2358 3.2437 – 3.2624 3.0160 3.0850 3.0845 3.0844 3.1005 3.0945 3.1009 2.8463a 2.8863a 10.6772 10.6387 10.6762 – 10.2370 10.2908 10.2178
0.91231 0.92097 0.91015 0.90357 – 0.8867 – 1.9129a 1.9652a 1.9451a 1.5265a 1.9312a 1.8505a 0.92671 0.93062 0.91253 0.84479 0.8896 0.8010 1.0273a 1.8232a 0.9110 0.90387 0.80674 – 1.8320a 2.1064a 0.89240 0.89617 0.90215 0.87393 – – 0.8647 1.9504a 0.90937 0.90648 0.90359 0.88868 – 0.8913 2.0720a 1.9755a 0.80326 0.84935 0.75042 – 1.7450a 1.2050a 1.8909a
0.13777 0.13786 0.13775 0.13833 – 0.1385 0.1382 0.1273 0.1261 – – 0.1264 0.1277 0.13812 0.13776 0.13832 0.13867 – – 0.1378 0.1282 0.13632 0.13671 0.13743 0.1366 0.1257 0.1243 0.13844 0.13844 0.13830 0.13888 0.12734 – 0.1387 0.1269 0.13836 0.13818 0.13832 0.13858 0.12604 0.1386 – 0.1265 0.14039 0.13934 0.14157 – 0.1289 0.1363 –
Case 2
Case 3
Case 4
Case 5
Case 6
a
Infeasible solution, constraint on load bus voltage is violated.
Fig. 8. Comparative convergence of case 3 for IEEE 30-bus system.
Fig. 7. Comparative convergence of case 2 for IEEE 30-bus system.
92
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 8 Comparison of 3-CH techniques with previous studies for multi-objective cases of IEEE 30-bus system. Case #
Algorithm
Fitness
Fuel cost ($/h)
Emission (t/h)
𝑃loss (MW)
VD (p.u.)
L-index (max)
Case 7
ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) MFO (Mohamed et al., 2017) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) MFO (Mohamed et al., 2017) BSA (Chaib et al., 2016) ICBO (Bouchekara et al., 2016b) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) FPA (Mohamed et al., 2017) BSA (Chaib et al., 2016) ICBO (Bouchekara et al., 2016b) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) MFO (Mohamed et al., 2017)
1040.151 1040.125 1040.134 1040.808 1041.671 813.1742 813.1956 813.1959 814.1545 814.3541 814.8994 813.5378 814.1708 814.1649 814.1841 814.9378 814.9067 812.9240a 811.8477a 964.1331 964.1254 964.1234 965.2905 965.8077
858.867 859.1458 858.9319 859.1915 858.5812 803.7198 803.4241 803.4196 803.3125 803.7911 803.4294 803.3978 800.4321 800.4203 800.4365 801.2248 801.1487 800.3340a 799.3277a 830.1156 830.1366 830.2123 830.639 830.9135
0.22902 0.2289 0.22895 0.22899 0.22947 0.36384 0.36424 0.36324 0.36344 0.36355 0.3546 – 0.36585 0.36592 0.36517 0.36106 0.3718 0.3514 – 0.25293 0.25313 0.25294 0.25258 0.25231
4.5321 4.5245 4.5301 4.5404 4.5772 9.8414 9.7807 9.7573 9.7206 9.8685 9.3751 9.7453 9.0043 8.9985 8.9838 8.9761 9.3174 8.4904 8.6465 5.5894 5.5887 5.5857 5.6219 5.5971
0.93028 0.92731 0.92626 0.92852 0.89944 0.09454 0.09772 0.09776 0.10842 0.10563 0.1147 0.1014 0.91244 0.93846 0.93706 0.92655 0.87563 1.9855a 1.9961a 0.29738 0.29653 0.29615 0.29385 0.33164
0.13796 0.13785 0.13781 0.13814 0.13806 0.14888 0.14893 0.14893 0.14783 0.14906 0.14840 0.14900 0.13739 0.13745 0.13748 0.13713 0.13758 0.12590 0.12520 0.14748 0.14756 0.14756 0.14802 0.14556
Case 8
Case 9
Case 10
a
Infeasible solution, constraint on load bus voltage is violated.
Fig. 11. Comparative convergence of case 6 for IEEE 30-bus system.
Fig. 9. Comparative convergence of case 4 for IEEE 30-bus system.
Fig. 12. Convergence of case 7 (SF-DE) for IEEE 30-bus system. Fig. 10. Comparative convergence of case 5 for IEEE 30-bus system.
5.3. System 2: IEEE 57-bus test system value is included in the diagram. As an obvious fact, fitness function Settings of control and state variables alongwith their ranges for the best fitness values of study cases for 57-bus test system are listed in Table 9. In Table 10, the results of all 3-CH techniques (ECHT, SF and SP) are compared with most recent studies available for same objective functions. As it turns out, SP-DE leads to best fitness values in all the
gradually reduces in all these cases, though the individual objectives show irregularities during the process of convergence. The fact is due to non-linear relationship of the individual objective functions with control parameters of the network. 93
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 9 Simulation results of best solutions for system 2: IEEE 57-bus system. Parameters Algorithm 𝑃G2 (MW) 𝑃G3 (MW) 𝑃G6 (MW) 𝑃G8 (MW) 𝑃G9 (MW) 𝑃G12 (MW) 𝑉1 (p.u.) 𝑉2 (p.u.) 𝑉3 (p.u.) 𝑉6 (p.u.) 𝑉8 (p.u.) 𝑉9 (p.u.) 𝑉12 (p.u.) 𝑄c18 (MVAr) 𝑄c25 (MVAr) 𝑄c53 (MVAr) 𝑇19 (p.u.) 𝑇20 (p.u.) 𝑇31 (p.u.) 𝑇35 (p.u.) 𝑇36 (p.u.) 𝑇37 (p.u.) 𝑇41 (p.u.)
Min 30 40 30 100 30 100 0.95 0.95 0.95 0.95 0.95 0.95 0.95 0 0 0 0.90 0.90 0.90 0.90 0.90 0.90 0.90
Max
Case 11
Case 12
Case 13
Case 14
Parameters
Min
Max
Case 11
Case 12
Case 13
Case 14
100 140 100 550 100 410 1.10 1.10 1.10 1.10 1.10 1.10 1.10 20 20 20 1.10 1.10 1.10 1.10 1.10 1.10 1.10
SP-DE 89.759 44.564 71.5502 459.8236 98.2124 359.4297 1.0645 1.062 1.0548 1.0598 1.0754 1.0488 1.0481 7.556 13.2503 11.4127 1.046 0.9656 1.0192 1.0023 1.0356 1.0294 0.9955
SP-DE 89.1352 45.135 71.6172 460.5427 96.338 360.8131 1.0328 1.0311 1.0262 1.0414 1.0624 1.0277 1.0179 4.9653 14.9616 14.9965 0.9982 0.9796 0.9918 1.0358 0.9643 1.0239 1.0151
SF-DE 90.2581 44.907 71.1695 461.2539 94.1207 360.8523 1.062 1.0598 1.0528 1.0591 1.0752 1.0477 1.0473 7.5494 13.3646 12.0559 0.9971 0.9803 1.0104 1.0088 1.0149 1.0315 0.9964
SP-DE 68.1397 92.9072 34.0623 294.8929 95.3897 373.4152 1.0086 1.0065 1.0132 1.0021 1.0176 1.0082 1.0392 4.0015 19.7885 19.9272 0.923 1.0505 0.9712 1.0577 1.0624 1.0043 0.9926
𝑇46 (p.u.) 𝑇54 (p.u.) 𝑇58 (p.u.) 𝑇59 (p.u.) 𝑇65 (p.u.) 𝑇66 (p.u.) 𝑇71 (p.u.) 𝑇73 (p.u.) 𝑇76 (p.u.) 𝑇80 (p.u.) Fuel cost ($/h) Emission (t/hr) 𝑃loss (MW) VD (p.u.) L-index (max) 𝑃G1 (MW) 𝑄G1 (MVAr) 𝑄G2 (MVAr) 𝑄G3 (MVAr) 𝑄G6 (MVAr) 𝑄G8 (MVAr) 𝑄G9 (MVAr) 𝑄12 (MVAr) CPU time (s)
0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.90
1.10 1.10 1.10 1.10 1.10 1.10 1.10 1.10 1.10 1.10
0 −140 −17 −10 −8 −140 −3 −150
576 200 50 60 25 200 9 155
0.9543 0.9134 0.9783 0.9696 0.9754 0.9421 0.9743 0.9952 0.9595 0.9987 41667.82 1.350 14.9090 1.54367 0.28123 142.37 47.7018 48.6597 35.7354 −6.5273 54.45 8.7382 53.0263 219.9
0.9403 0.9015 0.966 0.9667 0.982 0.9349 0.9666 1.0007 0.9411 1.0075 41697.50 1.3550 15.5897 0.77253 0.29228 142.8084 42.0835 49.8259 33.8151 −4.9521 75.3832 8.8451 44.529 203.6
0.9638 0.9127 0.9772 0.9651 0.975 0.9397 0.9734 0.9911 0.9637 0.9921 41667.53 1.3576 14.8963 1.61174 0.28022 143.1347 44.3579 49.5019 33.7255 −6.6773 56.1475 8.2295 54.8241 214.4
0.9236 0.9002 0.9302 0.9866 1.0127 0.9002 0.9628 1.0016 0.9059 0.9895 45549.49 1.2898 18.4275 0.59267 0.30052 310.4204 −29.4561 41.862 58.879 −4.8572 32.9393 8.4604 149.6131 212.5
Fig. 13. Convergence of case 8 (ECHT-DE) for IEEE 30-bus system. Fig. 15. IEEE 57-bus system – voltage profiles of load buses for best solutions of case 11 to case 14.
value of cumulative voltage deviation (VD in Eq. (23)) would be 3 p.u. (i.e. 50 × 0.06 p.u.). Reported VD value in one case is found to be higher than 3 p.u. and the same is marked with a footnote in comparison Table 10. The reference paper pertaining to the case study adopted most common static penalty function approach for constraint handling. Case 11 of minimizing basic objective of fuel cost leads to a value of 41667.82 $/h by SP-DE, the lowest when compared with other valid results of recent studies as depicted in Table 10. Further, network loss is also among the best (14.909 MW) with the recommended settings of control variables by SP-DE algorithm. Study in case 12 of simultaneous minimization of fuel cost and voltage deviation (𝑉 𝐷) by SP-DE achieves fitness value close to the value obtained by DSA (Abaci and Yamacli, 2016) method. Significant improvement in fuel cost is observed (by SFDE) in multi-objective optimization of case 13 where both cost and 𝐿index (max) are minimized. Case 14 is single objective optimization of minimizing voltage deviation. The value reported by SP-DE is much better than the comparable study by APFPA (Mahdad and Srairi, 2016), though with the sacrifice of fuel cost which increases from previous studies following the recommended settings of control variables. Roy and Paul (2015) and Shaheen et al. (2017) have used higher upper limits for shunt VAR compensators (up to 30 MVAr); hence straight comparison with present study is not valid. In general, fitness values
Fig. 14. Convergence of case 9 (SF-DE) for IEEE 30-bus system.
cases except case 13. Similar to system 1 of 30-buses, constraints are all diligently satisfied by the constraint handling methods. Ranges for generator reactive power are derived from Zimmerman et al. (0000) and found to be quite narrow for some generators. The CH techniques duly satisfy the limits as listed values of generator reactive power suggest. Further, the permissible range of load bus voltages for 57-bus system is between [0.94 p.u.−1.06 p.u.]. Theoretically, if all 50-load buses of the system operate at the limits, the maximum possible absolute 94
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Fig. 19. Comparative convergence of case 14 for IEEE 57-bus system.
Fig. 16. Comparative convergence of case 11 for IEEE 57-bus system.
Fig. 17. Convergence of case 12 (SP-DE) for IEEE 57-bus system. Fig. 20. IEEE 118-bus system – voltage profiles of load buses for study cases with SP-DE.
Fig. 18. Convergence of case 13 (SF-DE) for IEEE 57-bus system. Fig. 21. Convergences of case 15 and case 16 (SP-DE) for IEEE 118-bus system.
of the objective functions achieved by SP-DE and SF-DE in all the cases are better than most other reported values in past literatures. Curves in Fig. 15 describe the voltage profiles of load buses for the best solutions of case studies performed for 57-bus system. Clearly all voltages are within allowable limits. However, in a couple of cases many bus voltages are found to be very close to the upper limit. The fact reiterates the need of extra caution in satisfying voltage constraint so that no system bus experiences overvoltage. Convergence curves applying 3-CH techniques are superimposed for case 11 in Fig. 16 which shows reasonably fast convergence to the fitness value. SP-DE appears to have converged faster than the other methods. All the curves are shown starting after a certain number of function
evaluations. It is worthwhile to mention that, during initial phase of the search process, the algorithm (CH technique) looks for feasible solutions. The actual convergence to optimal solution starts when the search process enters the feasible region. For 30-bus system, the CH techniques easily find the feasible region. But possibly due to narrow limits of generator reactive power in 57-bus system, the algorithm takes few hundreds (or sometimes few thousands) of evaluations to find the feasible zone. The convergence curves do not indicate the initial phase of search process when the fitness progression is erratic. Now, for twoobjective cases, convergence curves of only best performing method are drawn in case 12 and case 13 to maintain clarity in diagrams. 95
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table 10 Comparison of 3-CH techniques with previous studies for various cases of IEEE 57-bus system. Case #
Algorithm
Fitness
Fuel cost ($/h)
Emission (t/h)
𝑃loss (MW)
VD (p.u.)
L-index (max)
Case 11
ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) DSA (Abaci and Yamacli, 2016) ICBO (Bouchekara et al., 2016b) ARCBBO (Kumar and Premalatha, 2015) APFPA (Mahdad and Srairi, 2016) LTLBO (Ghasemi et al., 2015) MICA-TLA (Ghasemi et al., 2014) DE (Shaheen et al., 2017) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) DSA (Abaci and Yamacli, 2016) MFO (Mohamed et al., 2017) MICA-TLA (Ghasemi et al., 2014) ECHT-DE SF-DE SP-DE MSA (Mohamed et al., 2017) DSA (Abaci and Yamacli, 2016) MFO (Mohamed et al., 2017) ECHT-DE SF-DE SP-DE APFPA (Mahdad and Srairi, 2016) KHA (Roy and Paul, 2015) DE (Shaheen et al., 2017)
41670.56 41667.85 41667.82 41673.72 41686.82 41697.33 41686 41628.75a 41679.55 41675.05 41682 41776.48 41775.09 41774.75 41782.80 41775.60 41786.66 42013.08 41699.25 41695.55 41696.54 41703.48 41785.05 41707.66 0.60416 0.59584 0.59267 0.8909 0.5810b 0.5839b
41670.56 41667.85 41667.82 41673.72 41686.82 41697.33 41686 41628.75a 41679.55 41675.05 41682 41694.82 41697.52 41697.50 41714.98 41699.40 41718.87 41959.18 41671.09 41667.53 41668.45 41675.99 41761.22 41680.19 46813.22 45246.02 45549.49 43485.93 42006.44 –
1.36230 1.35816 1.350 1.9526 – – – – – – – 1.3597 1.35769 1.3550 1.9551 – 2.0149 – 1.3609 1.3576 1.35524 1.9188 – 1.9192 1.3379 1.23453 1.2898 – – –
14.9479 14.8864 14.9090 15.0526 – 15.5470 15.3769 14.0470 15.1589 15.0149 – 15.5806 15.5616 15.5897 15.9214 – 16.2189 19.909 15.0275 14.8963 15.012 15.0026 – 15.1026 19.0821 18.4697 18.4275 12.1513 – –
1.50319 1.64209 1.54367 1.5508 1.0833 1.3173 – 3.5571a – 1.6161 – 0.81659 0.77572 0.77253 0.6782 0.7620 0.6780 0.5390 1.56188 1.61174 1.60803 1.7236 1.0573 1.7245 0.60416 0.59584 0.59267 0.8909 0.5810b 0.5839b
0.28886 0.27971 0.28123 0.28392 0.24353 0.27760 – – – – – 0.29198 0.29262 0.29228 0.29533 0.2471 0.29525 – 0.28152 0.28022 0.28092 0.27481 0.2383 0.27467 0.3008 0.30135 0.30052 – 0.2985 –
Case 12
Case 13
Case 14
a
Infeasible solution, constraint on load bus voltage is violated.
b
Higher limits for the shunt compensators have been considered.
As expected, unevenness in individual objectives, especially voltage deviation, along the optimization path is observed for multi-objective cases represented in Figs. 17 and 18. Single objective of minimizing voltage deviation requires higher number of function evaluations for all methods to converge as seen in Fig. 19. It might be because of highly non-linear relationship of bus voltages with numerous control variables of 57-bus system.
5.5. Additional study case, AC1 – handling discrete variables This section of the study considers cost minimization of IEEE-30 bus system considering discrete variables in the optimization process. In DE, discrete variables can easily be accounted for by ‘rounding off ‘operation. Here, both (shunt) capacitor banks and transformer taps are considered discrete variables. The capacitor can be switched in at discrete steps of 200 kVAr (i.e. 0.2 MVAr). Transformer tap is considered in p.u. within range 0.90 p.u. to 1.10 p.u. The possible transformer taps considered in this case study are [0.90–0.92–0.94-. . . −1.08–1.10] p.u. During the search process, the control variables on shunt capacitor ratings and transformer taps are simply rounded off (up to required decimal points) to the nearest multiples of 0.2 (MVAr) and 0.02 (p.u.), respectively. The optimization case on generation cost is run with all the three algorithms i.e. ECHT-DE, SF-DE and SP-DE. SF-DE algorithm achieves the minimum cost value of 800.4352 $/h, which is slightly higher than the best value obtained in same case with all continuous variables. The complete results of all algorithms with control variables and other calculated parameters are provided in Table A.7 in Appendix.
5.4. System 3: IEEE 118-bus test system To verify effectiveness of the proposed method, the large-scale IEEE 118-bus system is considered for study purpose. In general, performance of SP-DE is found to be superior for higher number of variables in constrained optimization problems. Hence, SP-DE algorithm is applied to the system to minimize fuel cost (case 15) and power loss (case 16) independently. A total of 300,000 function evaluations (with execution time of about 40 min) are performed for a study case with population size (𝑁𝑃 ) of 30. The details of calculation of these objective functions are mentioned in Section 3. The control variables and calculated parameters pertaining to the best fitness values obtained in 5-trial runs of each case are provided in Tables A.5 and A.6 in Appendix. Allowable ranges of generator real power (Zimmerman et al., 0000) are explicitly mentioned in the table. Shunt compensator ratings are selectable within 0 to 25 MVAr. Permissible ranges of generator bus voltages and transformer tap settings are within [0.95, 1.10] and [0.90, 1.10], respectively. The optimal value of fuel cost in case 15 is 135055.7 $/h, and that of real power loss in case 16 is 17.6946 MW. The system has a large number of control variables and constraints, and the number of critical constraints on generator reactive power and load bus accounts up to 118 nos. All the constraints are duly satisfied by the algorithm. Load bus voltage profiles for both the study cases are provided in Fig. 20. Convergence characteristics are shown in Fig. 21. Due to presence of huge number of constraints, the algorithm takes about 25 000 function evaluations to find first set of feasible solutions. The convergence to optimal solutions is observed in about 200 000 function evaluations.
6. Conclusion This paper discusses in detail the application and usefulness of three constraint handling techniques – SF, SP and ensemble of these two with DE being the basic search algorithm for optimal power flow (OPF) problem which consists of non-linear objectives and non-linear constraints. A comparative study among the three techniques points out the difficulty in guaranteeing superiority of a CH technique over others under various scenarios of a real-world problem. Ensemble method helps to an extent in achieving near best solution almost in all cases. However, the ensemble method does not always promise most optimal solution and fast convergence. Nonetheless, the importance of an efficient constraint handling technique is paramount. As substantiated in our study, without proper constraint handling method, and especially with static 96
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table A.1 Branch data used for 30-bus system. Branch no.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
Bus no.
R (p.u.)
From
To
1 1 2 3 2 2 4 5 6 6 6 6 9 9 4 12 12 12 12 14 16
2 3 4 4 5 6 6 7 7 8 9 10 11 10 12 13 14 15 16 15 17
0.0192 0.0452 0.0570 0.0132 0.0472 0.0581 0.0119 0.0460 0.0267 0.0120 0 0 0 0 0 0 0.1231 0.0662 0.0945 0.2210 0.0524
X (p.u.)
B/2 (p.u.)
0.0575 0.1652 0.1737 0.0379 0.1983 0.1763 0.0414 0.1160 0.0820 0.0420 0.2080 0.5560 0.2080 0.1100 0.2560 0.1400 0.2559 0.1304 0.1987 0.1997 0.1923
0.0264 0.0204 0.0184 0.0042 0.0209 0.0187 0.0045 0.0102 0.0085 0.0045 0 0 0 0 0 0 0 0 0 0 0
Rating (MVA)
Branch no.
130 130 65 130 130 65 90 70 130 32 65 32 65 65 65 65 32 32 32 16 16
Bus no.
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
From
To
15 18 19 10 10 10 10 21 15 22 23 24 25 25 28 27 27 29 8 6
18 19 20 20 17 21 22 22 23 24 24 25 26 27 27 29 30 30 28 28
R (p.u.)
X (p.u.)
B/2 (p.u.)
Rating (MVA)
0.1073 0.0639 0.0340 0.0936 0.0324 0.0348 0.0727 0.0116 0.1000 0.1150 0.1320 0.1885 0.2544 0.1093 0 0.2198 0.3202 0.2399 0.0636 0.0169
0.2185 0.1292 0.0680 0.2090 0.0845 0.0749 0.1499 0.0236 0.2020 0.1790 0.2700 0.3292 0.3800 0.2087 0.3960 0.4153 0.6027 0.4533 0.2000 0.0599
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0214 0.0065
16 16 32 32 32 32 32 32 16 16 16 16 16 16 65 16 16 16 32 32
Table A.2 Cost and emission coefficients of generators for IEEE 30-bus system (Chaib et al., 2016; Alsac and Stott, 1974). Generator
Bus
𝑎
𝑏
𝑐
𝑑
𝑒
𝛼
𝛽
𝛾
𝜔
μ
𝐺1 𝐺2 𝐺3 𝐺4 𝐺5 𝐺6
1 2 5 8 11 13
0 0 0 0 0 0
2 1.75 1 3.25 3 3
0.00375 0.0175 0.0625 0.00834 0.025 0.025
18 16 14 12 13 13.5
0.037 0.038 0.04 0.045 0.042 0.041
4.091 2.543 4.258 5.326 4.258 6.131
−5.554 −6.047 −5.094 −3.55 −5.094 −5.555
6.49 5.638 4.586 3.38 4.586 5.151
0.0002 0.0005 0.000001 0.002 0.000001 0.00001
2.857 3.333 8 2 8 6.667
Table A.3 Cost coefficients of generators 1 and 2 for IEEE 30-bus system for multi-fuel sources (Chaib et al., 2016). Generator
Bus
𝑎
𝐺1 𝐺2
1 2
55 40
𝑏
𝑐
𝑃𝐺𝑚𝑖𝑛
𝑃𝐺𝑚𝑎𝑥
𝑎
50 20
140 55
82.5 80
𝑖1
0.7 0.3
0.005 0.01
𝑏
𝑐
𝑖1
𝑃𝐺𝑚𝑖𝑛
𝑃𝐺𝑚𝑎𝑥
140 55
200 80
𝑖2
1.05 0.6
0.0075 0.02
𝑖2
Table A.4 Cost and emission coefficients of generators for IEEE 57-bus system (Chaib et al., 2016; Zimmerman et al., 0000). Generator
Bus
𝑎
𝑏
𝑐
𝑑
𝑒
𝛼
𝛽
𝛾
𝜔
μ
𝐺1 𝐺2 𝐺3 𝐺4 𝐺5 𝐺6 𝐺7
1 2 3 6 8 9 12
0 0 0 0 0 0 0
20 40 20 40 20 40 20
0.0775795 0.01 0.25 0.01 0.0222222 0.01 0.0322581
18 16 13.5 18 14 15 12
0.037 0.038 0.041 0.037 0.04 0.039 0.045
4.091 2.543 6.131 3.491 4.258 2.754 5.326
−5.554 −6.047 −5.555 −5.754 −5.094 −5.847 −3.555
6.49 5.638 5.151 6.39 4.586 5.238 3.38
0.0002 0.0005 0.00001 0.0003 0.000001 0.0004 0.002
0.286 0.333 0.667 0.266 0.8 0.288 0.2
help faster convergence and betterment of solutions especially for largescale power systems.
penalty function approach, limits on network parameters may often unknowingly be violated. Violating physical or security constraints of network components may compromise system safety, lead to excessive losses, malfunction and often failure of the component. So, operating the network within defined limits is the pre-requisite for secure and proper operation. Recommended settings by the CH techniques successfully lead to the desired condition of the network. Furthermore, the proposed algorithms in this study perform noticeably better than many other equivalent optimization methods in finding solutions of OPF. Reduction in hourly operation cost has been established almost in all the cases studied under the scope of this literature. In future, the authors propose integration of advanced variants of DE with CH techniques that might
Acknowledgment This project is funded by the National Research Foundation Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) program. Appendix A See Table A.1–A.7. 97
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Table A.5 Simulation results of SP-DE algorithm for Case 15 of IEEE 118-bus system. Control variables
Range (MW)
Values (MW)
Control variables
Range (MW)
Values (MW)
Control variables
Values (p.u.)
Control variables
Values (p.u.)
Control variables
Values
𝑃G1 (MW)
30–100
30.0317
𝑃G65
147.3–491
289.6489
𝑉G1 (MW)
0.9871
𝑉G65
1.0623
𝑄C5 (MVAr)
14.0629
𝑃G4
30–100
30.0143
𝑃G66
147.6–492
289.1504
𝑉G4
1.0153
𝑉G66
1.0593
𝑄C34 (MVAr)
0.8808
𝑃G6
30–100
30.0122
𝑃G70
30–100
30.0116
𝑉G6
1.008
𝑉G69
1.0389
𝑄C37 (MVAr)
0.8844
𝑃G8
30–100
30.0052
𝑃G72
30–100
30.0013
𝑉G8
1.0388
𝑉G70
1.0195
𝑄C44 (MVAr)
5.768
𝑃G10
165–550
317.0191
𝑃G73
30–100
30.0033
𝑉G10
1.0494
𝑉G72
1.0191
𝑄C45 (MVAr)
21.5888
𝑃G12
55.5–185
66.8229
𝑃G74
30–100
30.0088
𝑉G12
1.002
𝑉G73
1.0234
𝑄C46 (MVAr)
10.9322
𝑃G15
30–100
30.0162
𝑃G76
30–100
30.0074
𝑉G15
0.9988
𝑉G74
1.0058
𝑄C48 (MVAr)
4.6786
𝑃G18
30–100
30.0063
𝑃G77
30–100
30.0141
𝑉G18
0.9989
𝑉G76
0.9868
𝑄C74 (MVAr)
24.2029
𝑃G19
30–100
30.0495
𝑃G80
173.1–577
350.9989
𝑉G19
0.9986
𝑉G77
1.013
𝑄C79 (MVAr)
23.8787
𝑃G24
30–100
30.0111
𝑃G85
30–100
30.0087
𝑉G24
1.015
𝑉G80
1.0218
𝑄C82 (MVAr)
23.5807
𝑃G25
96–320
152.1726
𝑃G87
31.2–104
31.2015
𝑉G25
1.0298
𝑉G85
1.0242
𝑄C83 (MVAr)
20.4897
𝑃G26
124.2–414
220.8106
𝑃G89
212.1–707
379.9452
𝑉G26
1.0744
𝑉G87
1.0432
𝑄C105 (MVAr)
13.4731
𝑃G27
30–100
30.0364
𝑃G90
30–100
30.0443
𝑉G27
1.0045
𝑉G89
1.0274
𝑄C107 (MVAr)
1.936
𝑃G31
32.1–107
32.1004
𝑃G91
30–100
30.021
𝑉G31
0.9992
𝑉G90
1.0062
𝑄C110 (MVAr)
18.1676
𝑃G32
30–100
30.0143
𝑃G92
30–100
30.0162
𝑉G32
1.0045
𝑉G91
1.0074
𝑇8 (p.u.)
1.0148
𝑃G34
30–100
30.0024
𝑃G99
30–100
30.0027
𝑉G34
1.0141
𝑉G92
1.0153
𝑇32 (p.u.)
1.0978
𝑃G36
30–100
30.0132
𝑃G100
105.6–352
177.1013
𝑉G36
1.0105
𝑉G99
1.0182
𝑇36 (p.u.)
1.0348
𝑃G40
30–100
30.0216
𝑃G103
42–140
42.0053
𝑉G40
1.0001
𝑉G100
1.0187
𝑇51 (p.u.)
1.0107
𝑃G42
30–100
30.038
𝑃G104
30–100
30.0088
𝑉G42
1.0081
𝑉G103
1.0146
𝑇93 (p.u.)
0.9946
𝑃G46
35.7–119
35.7003
𝑃G105
30–100
30.0022
𝑉G46
1.0316
𝑉G104
1.0067
𝑇95 (p.u.)
1.0095
𝑃G49
91.2–304
162.3848
𝑃G107
30–100
30.013
𝑉G49
1.0429
𝑉G105
1.0063
𝑇102 (p.u.)
0.9771
𝑃G54
44.4–148
44.6599
𝑃G110
30–100
30.0043
𝑉G54
1.0217
𝑉G107
0.9993
𝑇107 (p.u.)
0.9715
𝑃G55
30–100
30.0485
𝑃G111
40.8–136
40.8014
𝑉G55
1.0215
𝑉G110
1.0147
𝑇127 (p.u.)
1.0214
𝑃G56
30–100
30.0079
𝑃G112
30–100
30.0166
𝑉G56
1.0215
𝑉G111
1.0247
Fuel cost ($/h)
135055.7
𝑃G59
76.5–255
125.3306
𝑃G113
30–100
30.0223
𝑉G59
1.0424
𝑉G112
1.0046
𝑃loss (MW)
60.9596
𝑃G61
78–260
124.1197
𝑃G116
30–100
30.0052
𝑉G61
1.0496
𝑉G113
1.0057
VD (p.u.)
1.0715
𝑃G62
30–100
30.0168
𝑉G62
1.0462
𝑉G116
1.0592
𝑃G69 (MW) [0 – 805.2]
370.4283
Values (MW)
Control variables
Values (p.u.)
Control variables
Values (p.u.)
Control variables
Values
Table A.6 Simulation results of SP-DE algorithm for Case 16 of IEEE 118-bus system. Control variables
Range (MW)
Values (MW)
Control variables
Range (MW)
𝑃G1 (MW)
30–100
67.1484
𝑃G65
147.3–491
147.3227
𝑉G1 (MW)
0.9888
𝑉G65
1.0134
𝑄C5 (MVAr)
9.3906
𝑃G4
30–100
30.0136
𝑃G66
147.6–492
147.606
𝑉G4
1.005
𝑉G66
1.047
𝑄C34 (MVAr)
0.322
𝑃G6
30–100
30.95
𝑃G70
30–100
30.1083
𝑉G6
1.0009
𝑉G69
1.0177
𝑄C37 (MVAr)
0.0376
𝑃G8
30–100
30.0116
𝑃G72
30–100
30.0151
𝑉G8
1.029
𝑉G70
1.0149
𝑄C44 (MVAr)
3.1065
𝑃G10
165–550
165.0121
𝑃G73
30–100
30.029
𝑉G10
1.04
𝑉G72
1.0215
𝑄C45 (MVAr)
21.2411
𝑃G12
55.5–185
136.4669
𝑃G74
30–100
95.4903
𝑉G12
1
𝑉G73
1.02
𝑄C46 (MVAr)
5.8559
𝑃G15
30–100
84.631
𝑃G76
30–100
99.8349
𝑉G15
1.0008
𝑉G74
1.0127
𝑄C48 (MVAr)
8.9501
𝑃G18
30–100
32.5503
𝑃G77
30–100
99.9777
𝑉G18
0.9997
𝑉G76
0.9992
𝑄C74 (MVAr)
23.8294
𝑃G19
30–100
58.8217
𝑃G80
173.1–577
289.4337
𝑉G19
0.9992
𝑉G77
1.0079
𝑄C79 (MVAr)
24.9914
𝑃G24
30–100
30.0197
𝑃G85
30–100
31.0829
𝑉G24
1.0164
𝑉G80
1.0131
𝑄C82 (MVAr)
24.8441
𝑃G25
96–320
96.0079
𝑃G87
31.2–104
31.2094
𝑉G25
1.0239
𝑉G85
1.0124
𝑄C83 (MVAr)
14.0215
𝑃G26
124.2–414
124.2228
𝑃G89
212.1–707
212.105
𝑉G26
1.0443
𝑉G87
1.0332
𝑄C105 (MVAr)
11.5006
𝑃G27
30–100
49.1688
𝑃G90
30–100
99.7784
𝑉G27
1.0062
𝑉G89
1.0181
𝑄C107 (MVAr)
20.4597
𝑃G31
32.1–107
61.5931
𝑃G91
30–100
30.0022
𝑉G31
1.0048
𝑉G90
1.0096
𝑄C110 (MVAr)
17.3712
𝑃G32
30–100
39.6414
𝑃G92
30–100
30.1177
𝑉G32
1.0058
𝑉G91
1.0128
𝑇8 (p.u.)
1.0224
𝑃G34
30–100
63.2896
𝑃G99
30–100
38.6397
𝑉G34
1.0055
𝑉G92
1.0093
𝑇32 (p.u.)
1.0727
𝑃G36
30–100
53.3235
𝑃G100
105.6–352
105.6473
𝑉G36
1.0015
𝑉G99
1.0114
𝑇36 (p.u.)
1.0129
𝑃G40
30–100
99.9392
𝑃G103
42–140
42.1357
𝑉G40
1.0087
𝑉G100
1.0116
𝑇51 (p.u.)
1.0044
𝑃G42
30–100
99.9191
𝑃G104
30–100
30.8949
𝑉G42
1.0175
𝑉G103
1.0138
𝑇93 (p.u.)
0.9599
𝑃G46
35.7–119
82.1758
𝑃G105
30–100
52.9999
𝑉G46
1.0371
𝑉G104
1.0117
𝑇95 (p.u.)
0.9583
𝑃G49
91.2–304
141.0635
𝑃G107
30–100
57.6295
𝑉G49
1.0395
𝑉G105
1.0123
𝑇102 (p.u.)
0.9431
𝑃G54
44.4–148
147.8306
𝑃G110
30–100
30.1049
𝑉G54
1.0519
𝑉G107
1.012
𝑇107 (p.u.)
0.9677 0.9847
𝑃G55
30–100
79.5092
𝑃G111
40.8–136
40.8554
𝑉G55
1.052
𝑉G110
1.017
𝑇127 (p.u.)
𝑃G56
30–100
99.8073
𝑃G112
30–100
51.6249
𝑉G56
1.0509
𝑉G111
1.0247
Fuel cost ($/h)
155724.9
𝑃G59
76.5–255
251.0927
𝑃G113
30–100
30.0672
𝑉G59
1.05
𝑉G112
1.014
𝑃loss (MW)
17.6946
𝑃G116
30–100
75.151
𝑃G61
78–260
78.2118
𝑃G62
30–100
64.477
𝑉G61
1.0493
𝑉G113
1.0075
VD (p.u.)
0.8663
𝑉G62
1.0473
𝑉G116
1.0101
𝑃G69 (MW) [0 – 805.2]
2.9323
98
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100 Table A.7 Simulation results of study case AC1: handling discrete variables. Parameters
Range
ECHT-DE
SF-DE
SP-DE
𝑃G2 (MW) 𝑃G5 (MW) 𝑃G8 (MW) 𝑃G11 (MW) 𝑃G13 (MW) 𝑉1 (p.u.) 𝑉2 (p.u.) 𝑉5 (p.u.) 𝑉8 (p.u.) 𝑉11 (p.u.) 𝑉13 (p.u.) 𝑄c10 (MVAr) 𝑄c12 (MVAr) 𝑄c15 (MVAr) 𝑄c17 (MVAr) 𝑄c20 (MVAr) 𝑄c21 (MVAr) 𝑄c23 (MVAr) 𝑄c24 (MVAr) 𝑄c29 (MVAr) 𝑇11 (p.u.) 𝑇12 (p.u.) 𝑇15 (p.u.) 𝑇36 (p.u.) Fuel cost ($/h) Emission (t/h) 𝑃loss (MW) VD (p.u.) L-index (max) 𝑃G1 (MW)
20 – 80 15 – 50 10 – 35 10 – 30 12 – 40 0.95 – 1.10 0.95 – 1.10 0.95 – 1.10 0.95 – 1.10 0.95 – 1.10 0.95 – 1.10 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.0 – 5.0 0.90 – 1.10 0.90 – 1.10 0.90 – 1.10 0.90 – 1.10
48.4217 21.5846 21.3143 11.9471 12.0404 1.0821 1.0631 1.032 1.0357 1.0676 1.0533 3.2 1.2 4.4 5 3.4 5 3 5 2.8 1.04 0.92 0.98 0.98 800.4430 0.36601 8.9896 0.8689 0.1392 177.0814
48.5786 21.4729 21.0862 11.91 12.0122 1.0817 1.0627 1.0314 1.0363 1.0767 1.0532 1 1.4 4.4 5 3.6 5 2.8 5 2.8 1.04 0.92 0.98 0.98 800.4352 0.36684 9.0179 0.87519 0.13906 177.3581
48.36621 21.34167 22.08471 11.80553 12.03479 1.083016 1.064776 1.0343 1.039547 1.08065 1.037134 4.2 3.2 5 4.4 3.2 4.6 3.2 4.8 3 1.06 0.92 0.96 0.98 800.4478 0.36510 8.9648 0.88987 0.13860 176.7319
50 – 200
References
Deb, K., 2000. An efficient constraint handling method for genetic algorithms. Comput. Methods Appl. Mech. Eng. 186 (2), 311–338. Derrac, J., García, S., Molina, D., Herrera, F., 2011. A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms. Swarm Evol. Comput. 1 (1), 3–18. Farmani, R., Wright, J.A., 2003. Self-adaptive fitness formulation for constrained optimization. IEEE Trans. Evol. Comput. 7 (5), 445–455. Ghasemi, M., Ghavidel, S., Gitizadeh, M., Akbari, E., 2015. An improved teaching– learning-based optimization algorithm using Lévy mutation strategy for non-smooth optimal power flow. Int. J. Electr. Power Energy Syst. 65, 375–384. Ghasemi, M., Ghavidel, S., Rahmani, S., Roosta, A., Falah, H., 2014. A novel hybrid algorithm of imperialist competitive algorithm and teaching learning algorithm for optimal power flow problem with non-smooth cost functions. Eng. Appl. Artif. Intell. 29, 54–69. Kessel, P., Glavitsch, H., 1986. Estimating the voltage stability of a power system. IEEE Trans. Power Deliv. 1 (3), 346–354. Kumar, A.R., Premalatha, L., 2015. Optimal power flow for a deregulated power system using adaptive real coded biogeography-based optimization. Int. J. Electr. Power Energy Syst. 73, 393–399. Li, C., Yang, S., Nguyen, T.T., Yu, E.L., Yao, X., Jin, Y., Beyer, H.G., Suganthan, P.N., (2008). Benchmark generator for CEC 2009 competition on dynamic optimization. Liang, J.J., Qu, B.Y., Suganthan, P.N., (2013). Problem definitions and evaluation criteria for the CEC 2014 special session and competition on single objective real-parameter numerical optimization. Computational Intelligence Laboratory, Zhengzhou University, Zhengzhou China and Technical Report, Nanyang Technological University, Singapore. Liang, J.J., Runarsson, T.P., Mezura-Montes, E., Clerc, M., Suganthan, P.N., Coello, C.C., Deb, K., 2006. Problem definitions and evaluation criteria for the CEC 2006 special session on constrained real-parameter optimization. J. Appl. Mech. 41 (8). Mahdad, B., Srairi, K., 2016. Security constrained optimal power flow solution using new adaptive partitioning flower pollination algorithm. Appl. Soft Comput. 46, 501–522. Mallick, S., Acharjee, P., Ghoshal, S.P., Thakur, S.S., 2013. Determination of maximum load margin using fuzzy logic. Int. J. Electr. Power Energy Syst. 52, 231–246. Mallipeddi, R., Jeyadevi, S., Suganthan, P.N., Baskar, S., 2012. Efficient constraint handling for optimal reactive power dispatch problems. Swarm Evol. Comput. 5, 28– 36. Mallipeddi, R., Suganthan, P.N., 2010. Ensemble of constraint handling techniques. IEEE Trans. Evol. Comput. 14 (4), 561–579. Mohamed, A.A.A., Mohamed, Y.S., El-Gaafary, A.A., Hemeida, A.M., 2017. Optimal power flow using moth swarm algorithm. Electr. Power Syst. Res. 142, 190–206. Naderi, E., Narimani, H., Fathi, M., Narimani, M.R., 2017. A novel fuzzy adaptive configuration of particle swarm optimization to solve large-scale optimal reactive power dispatch. Appl. Soft Comput. 53, 441–456.
Abaci, K., Yamacli, V., 2016. Differential search algorithm for solving multi-objective optimal power flow problem. Int. J. Electr. Power Energy Syst. 79, 1–10. Acharjee, P., 2012. Identification of maximum loadability limit and weak buses using security constraint genetic algorithm. Int. J. Electr. Power Energy Syst. 36 (1), 40–50. Acharjee, P., Mallick, S., Thakur, S.S., Ghoshal, S.P., 2011. Detection of maximum loadability limits and weak buses using Chaotic PSO considering security constraints. Chaos Solut. Fractals 44 (8), 600–612. Alsac, O., Stott, B., 1974. Optimal load flow with steady-state security. IEEE Trans. Power Apparatus Syst. (3), 745–751. Ayan, K., Kılıç, U., Baraklı, B., 2015. Chaotic artificial bee colony algorithm based solution of security and transient stability constrained optimal power flow. Int. J. Electr. Power Energy Syst. 64, 136–147. Bhattacharya, A., Chattopadhyay, P.K., 2011. Application of biogeography-based optimisation to solve different optimal power flow problems. IET Gener. Transm. Distribution 5 (1), 70–80. Bhowmik, A.R., Chakraborty, A.K., 2015. Solution of optimal power flow using nondominated sorting multi objective opposition based gravitational search algorithm. Int. J. Electr. Power Energy Syst. 64, 1237–1250. Biswas, P.P., Suganthan, P.N., Amaratunga, G.A., 2017. Optimal power flow solutions incorporating stochastic wind and solar power. Energy Convers. Manag. 148, 1194– 1207. Bouchekara, H.R.E.H., Chaib, A.E., Abido, M.A., 2016a. Multiobjective optimal power flow using a fuzzy based grenade explosion method. Energy Syst. 7 (4), 699–721. Bouchekara, H.R.E.H., Chaib, A.E., Abido, M.A., El-Sehiemy, R.A., 2016b. Optimal power flow using an improved colliding bodies optimization algorithm. Appl. Soft Comput. 42, 119–131. Cai, H.R., Chung, C.Y., Wong, K.P., 2008. Application of differential evolution algorithm for transient stability constrained optimal power flow. IEEE Trans. Power Syst. 23 (2), 719–728. Chaib, A.E., Bouchekara, H.R.E.H., Mehasni, R., Abido, M.A., 2016. Optimal power flow with emission and non-smooth cost functions using backtracking search optimization algorithm. Int. J. Electr. Powerand Energy Syst. 81, 64–77. Chen, Q., Liu, B., Zhang, Q., Liang, J.J., Suganthan, P.N., Qu, B.Y., (2014). Problem definition and evaluation criteria for CEC 2015 special session and competition on bound constrained single-objective computationally expensive numerical optimization. Computational Intelligence Laboratory, Zhengzhou University, China and Nanyang Technological University, Singapore, Tech. Rep. Coello, C.A.C., 2002. Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: A survey of the state of the art. Comput. Methods Appl. Mech. Eng. 191 (11), 1245–1287. Daryani, N., Hagh, M.T., Teimourzadeh, S., 2016. Adaptive group search optimization algorithm for multi-objective optimal power flow problem. Appl. Soft Comput. 38, 1012–1024. 99
P.P. Biswas et al.
Engineering Applications of Artificial Intelligence 68 (2018) 81–100
Niknam, T., Azizipanah-Abarghooee, R., Narimani, M.R., 2012. A new multi objective optimization approach based on TLBO for location of automatic voltage regulators in distribution systems. Eng. Appl. Artif. Intell. 25 (8), 1577–1588. Pandiarajan, K., Babulal, C.K., 2016. Fuzzy harmony search algorithm based optimal power flow for power system security enhancement. Int. J. Electr. Power Energy Syst. 78, 72–79. Pulluri, H., Naresh, R., Sharma, V., 2016. A solution network based on stud krill herd algorithm for optimal power flow problems. Soft. Comput. 1–18. Pulluri, H., Naresh, R., Sharma, V., 2017. Application of stud krill herd algorithm for solution of optimal power flow problems. Int. Trans. Electr. Energy Syst. 27 (6). Qin, A.K., Huang, V.L., Suganthan, P.N., 2009. Differential evolution algorithm with strategy adaptation for global numerical optimization. IEEE Trans. Evol. Comput. 13 (2), 398–417. Rabiee, A., Parniani, M., 2013. Voltage security constrained multi-period optimal reactive power flow using benders and optimality condition decompositions. IEEE Trans. Power Syst. 28 (2), 696–708. Reddy, S.S., Rathnam, C.S., 2016. Optimal power flow using glowworm swarm optimization. Int. J. Electr. Power Energy Syst. 80, 128–139. Roy, P.K., Paul, C., 2015. Optimal power flow using krill herd algorithm. Int. Trans. Electr. Energy Syst. 25 (8), 1397–1419.
Roy, R., Jadhav, H.T., 2015. Optimal power flow solution of power system incorporating stochastic wind power using Gbest guided artificial bee colony algorithm. Int. J. Electr. Power Energy Syst. 64, 562–578. Sedighizadeh, M., Sarvi, M., Naderi, E., 2011. Multi objective optimal power flow with FACTS devices using shuffled frog leaping algorithm. Int. Rev. Electr. Eng. 6 (4). Shaheen, A.M., El-Sehiemy, R.A., Farrag, S.M., 2016. Solving multi-objective optimal power flow problem via forced initialised differential evolution algorithm. IET Gener. Transm. Distribution 10 (7), 1634–1647. Shaheen, A.M., Farrag, S.M., El-Sehiemy, R.A., 2017. MOPF solution methodology. IET Gener. Transm. Distribution 11 (2), 570–581. Storn, R., Price, K., 1997. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optim. 11 (4), 341–359. Tessema, B., Yen, G.G., 2006. A self-adaptive penalty function based algorithm for constrained optimization. In: IEEE Congress on Evolutionary Computation, 2006. CEC 2006. IEEE, pp. 246–253. Wolpert, D.H., Macready, W.G., 1997. No free lunch theorems for optimization. IEEE Trans. Evol. Comput. 1 (1), 67–82. Zimmerman, R.D., Murillo-Sánchez, C.E., Thomas, R.J., Matpower. Available at: http: //www.pserc.cornell.edu/matpower.
100