# A Heuristic Approach for Optimal Planning and Operation of Distribution Systems

Authors: Khalid Mohammed Saffer Alzaidi; Oguz Bayat; Osman N. Uçan
Journal: Journal of Optimization
 (2018)

---

## Abstract

The efficient planning and operation of power distribution systems are becoming increasingly significant with the integration of renewable energy options into power distribution networks. Keeping voltage magnitudes within permissible ranges is vital; hence, control devices, such as tap changers, voltage regulators, and capacitors, are used in power distribution systems. This study presents an optimization model that is based on three heuristic approaches, namely, particle swarm optimization, imperialist competitive algorithm, and moth flame optimization, for solving the voltage deviation problem. Two different load profiles are used to test the three modified algorithms on IEEE 123- and IEEE 13-bus test systems. The proposed optimization model uses three different cases: Case 1, changing the tap positions of the regulators; Case 2, changing the capacitor sizes; and Case 3, integrating Cases 1 and 2 and changing the locations of the capacitors. The numerical results of the optimization model using the three heuristic algorithms are given for the two specified load profiles.


---

## Body

## 1. Introduction

Power systems have been evolving in the last two decades, exhibiting such changes as deregulation and the integration of renewables into the philosophical and operational mentalities. From the operational point of view, control means that involving the coordinated operation of tap changing transformers, such as capacitors, is required because loads are not constant over time and the outputs of renewable energy sources are intermittent. Voltage optimization (VO) is an effective technology that has been saving the industry millions of dollars in wasted electrical energy since the beginning of the new millennium [1]. High demand used to be managed by voltage reduction [2]. Another way of helping system operation is using capacitors to improve the power factor and the voltage profile and reduce power losses [3]. Furthermore, tap operations of voltage regulators are helpful in enhancing voltage profiles. Capacitors and voltage regulators are integrated, and improved voltage profiles are obtained. However, the life span of these devices is shortened by frequent operation because they are based on mechanical switch operations. New technological developments have made electronics-based voltage regulators and capacitors available [4], thereby bringing additional flexibility into the operation of smart grids.On the planning side, optimal capacitor locations are sought [4]. For instance, in an algorithm that depends on dynamic programming, fuzzy logic and genetic algorithm (GA) approaches are used for capacitor distribution in distribution feeders. Gravitational search algorithm was used for optimal capacitor placement in [5], whereas a teaching-learning-based optimization was used for the same aim in [6]. Capacitors can also be used to reduce the effects of harmonics in distribution systems; the harmony search algorithm was applied for this goal in [7]. Capacitor location and sizing problem have been solved by other heuristics, such as clonal selection algorithm [8], ant colony optimization algorithm [8], and PSO [9].Producing the best possible result with the available resources is always an objective in engineering problems. Optimization problems are generally solved using two approaches. The first is based on mathematical analysis, and the second is based on numerical calculations. Numerical optimization methods can be divided into derivative-based and non-derivative-based methods. If the derivatives of the encountered model are not easy to find or a mathematical function related to the model does not exist, then non-derivative-based methods are applied. These methods are generally inspired by nature. The most popular model is GA, which reflects the evolution process in nature [10]. Subsequently, methods inspired by the behaviors of birds and fish (particle swarm optimization [PSO]) [11], improvisation process of musicians (harmony search) [12], and the navigation approach of moths in nature, which is named transverse orientation (moth flame optimization [MFO]) [13], were developed.This work models the voltage optimization problem using three different heuristic algorithms, namely, imperialist competitive algorithm (ICA), particle swarm optimization (PSO), and moth flame optimization (MFO). Cases1 and 2 are applicable to operation, and Case 3 is applicable to planning in distribution systems.(i)
The first model changes and uses the tap positions of the voltage regulators and obtains the optimal voltage value for given load conditions of the distribution system.
(ii)
The second model uses only the capacitors and optimizes the sizes of these devices for given load conditions.
(iii)
The third model uses the voltage regulators and the capacitors and finds the optimal tap positions, capacitor sizes, and locations.MATLAB and a free power distribution system simulation tool OpenDSS [14, 15] are used in the simulations.The rest of the paper is organized as follows. Section2 proposes the voltage optimization models. Section 3 briefly explains ICA, PSO, MFO, and modified algorithm-based voltage deviation. Section 4 presents the experiments and the simulation results. Section 5 presents the conclusions.


## 2. Model

We model three different cases.Case 1.
This case considers tap changers for the voltage regulators to minimize voltage deviations. The optimization model is as follows:(1)
M
i
n
i
m
i
z
e
x
=
∑
i
=
1
N
V
i
-
1
2
,
V
i
=
f
T
a
p
i
(2)
Subject  to:
0.95
≤
V
i
≤
1.05
(3)
T
a
p
m
i
n
≤
T
a
p
i
≤
T
a
p
m
a
x
,where x denotes the fitness values (cost), N is the number of buses, V
i is the voltage magnitude of bus i, Tapi is the tap position of the regulator, and T
a
p
m
i
n and T
a
p
m
a
x represent the minimum and maximum positions that a tap in a regulator can take, respectively. These values are in the range of [
-
16,16
].Case 2.
This case considers changing the size of the capacitors, and the model is as follows:(4)
M
i
n
i
m
i
z
e
x
=
∑
i
=
1
N
V
i
-
1
2
,
V
i
=
f
C
a
p
i
(5)
Subject  to:
0.95
≤
V
i
≤
1.05
0
≤
C
a
p
i
≤
C
a
p
m
a
x
,where x represents the fitness values (cost), N is the number of buses, V
i is the voltage magnitude of bus i, C
a
p
i is the size of the bank capacitor, and C
a
p
m
a
x is the maximum size of the bank capacitor.Case 3.
This case integrates Cases1 and 2 and changes the locations of the capacitors. The mathematical model is as follows: (6)
M
i
n
i
m
i
z
e
x
=
∑
i
=
1
N
V
i
-
1
2
,
V
i
=
f
T
a
p
i
,
C
a
p
i
,
l
i
(7)
Subject  to:
0.95
≤
V
i
≤
1.05
T
a
p
m
i
n
≤
T
a
p
i
≤
T
a
p
m
a
x
0
≤
C
a
p
i
≤
C
a
p
m
a
x
2
≤
l
i
≤
L
m
a
x
,where x represents the fitness values (cost), N is the number of buses, V
i is the voltage magnitude of bus i, T
a
p
i is the tap position of the regulator, T
a
p
m
i
n and T
a
p
m
a
x represent the minimum and maximum positions that a tap in a regulator can take, respectively (these values are in the range of [
-
16,16
]), C
a
p
i is the size of the bank capacitor, C
a
p
m
a
x is the maximum size of the bank capacitor, l
i represents the location of the capacitors, and L
m
a
x represents the maximum bus location.


## 3. Heuristic Algorithms

### 3.1. Imperialist Competitive Algorithm (ICA)

#### 3.1.1. General Approach

ICA was recently developed in 2007 by Esmaeil Gargari and Caro Lucas for continuous optimization problems [16]. The working philosophy corresponds to other evolutionary algorithms and initially creates random solution candidates called countries. The cost function of each solution candidate shows the power of each country. Hence, populations are composed of either colonized or imperialist countries. According to random rules, a part of a population is selected as the imperialists or the powerful countries, and the remaining part of the population comprises the colonized. Figure 1 presents a flowchart of ICA [16].Figure 1
Flowchart of ICA [16].The method is conducted as follows:(i)
Form countries: theith country is formed as follows: (8)
c
o
u
n
t
r
y
i
=
P
1
,
P
2
,
P
3
,
…
,
P
D
N
,
where DN denotes the problem variables or dimensions. Initial random values forP
i should be within the lower and upper ranges for each variable.
(ii)
Find the powers of each country by evaluating the objective function of the optimization problem as follows:(9)
f
c
o
u
n
t
r
y
i
=
f
P
1
,
P
2
,
P
3
,
…
,
P
D
N
(iii)
Select the imperialist and colonized countries. The power of a country is inversely symmetrical to its cost. The division of colonies among imperialists and the normalized value of each imperialist is defined as follows:(10)
C
n
=
c
n
-
max
⁡
c
i
,
wherec
n is the cost of nth imperialist and C
n is the normalized value.
(iv)
Then, the colony countries move to the imperialist ones to start the optimization process. The DN country population is generated, andN
i
m
p represents the most powerful population, whose members are selected as imperialists (the sets of controller coefficients with smaller cost function in this problem). The remaining N
c
o
l countries are the colonies (the sets of controller coefficients with a high cost function in this problem), each of which is a part of one of the above-mentioned empires. In the attraction policy, the colonies move toward the imperialists along Mx units and are situated in a new position. Mx is a random variable with regular distribution and can be expressed as follows: (11)
M
x
~
U
0
,
β
·
d
s
,
whereβ is a constant number greater than 1 (
β
=
2
) and ds is the space between imperialist and colony. To search different points around the imperialist, we add a random amount θ of deviation to the direction of movement as follows:
θ
~
U
-
ץ
,
ץ, where ץ is a parameter to adjust the deviation value (
ץ
=
π
/
4
).
(v)
Calculate total power of an empire. It can be determined by the power of imperialist country plus percentage of power of its colonies as follows:
T.C.n = cost (imperialist ) + Ƹ · mean (cost (colonies of empire)), where T.C.n is the total cost of nth empire and Ƹ is a positive number that is considered to be less than 1 (Ƹ
=
0.02). The weakest colony of the weakest empire is picked out.There are some other hyperparameters used in the internal calculations of this algorithm; for example, the percent of search space size is a positive number (0.02), which enables the uniting process of two empires,α is a number in the interval of 0
1 and denotes the importance of mean minimum compares to the global minimum, and revolution rate is a positive number (0.3) representing the process in which the sociopolitical characteristics of a country change suddenly.


#### 3.1.2. ICA-Based Voltage Deviation Algorithm

The flowchart of the modified ICA algorithms is shown in Figure2, and the steps are as follows.Figure 2
Flowchart of modified ICA.Step 1.
Initialize the ICA parameters, namely, population sizeN
p
o
p, maximum iteration number M
a
x
i
t, number of imperialist countries N
i
m
p, and number of colony countries N
c
o
l. Set the voltage magnitude limits, and set the possible capacitor locations, capacitor size limits, and minimum and maximum tap positions depending on the case being solved.Step 2.
Randomly create the size and location of the capacitors and tap positions of the regulators, and form the initial country as follows:(12)
c
o
u
n
t
r
y
i
=
C
a
p
1
,
…
,
C
a
p
m
,
l
1
,
…
,
l
m
,
T
a
p
1
,
…
,
T
a
p
n
,where m and n represent the numbers of bank capacitors and voltage regulators, respectively.Step 3.
Run a load flow using the specified load profile and the solution candidates, perform a power flow, and calculate the fitness value of the test system depending on the case number as in (1), (4), and (6).Step 4.
Determine the imperialist and colonized countries depending on the fitness value as in (9) and (10).Step 5.
Update the size and location of the capacitors and the tap position of the regulators for all empires as in (11).Step 6.
Repeat Steps3–5 until the stopping condition is met.


### 3.2. Particle Swarm Optimization (PSO)

#### 3.2.1. General Approach

PSO was originally developed in 1995 by Kennedy and Eberhart and inspired by the social behavior of schooling fish and flocking birds [17]. The birds in a group are considered an individual in the PSO method. These particles can be flown through a search space. The location of a particle in the search problem represents one solution for the problem. A new and different solution is created when the individual moves to a new location in the search space. Each solution can be evaluated using an objective function that supplies a cost of the benefit of the solution. The direction and velocity of each individual can move along all dimensions of the search space and thus can change with all generation of movement. PSO is generally considered an evolutionary computation (EC) sample. Other EC examples include evolutionary strategies, genetic programming, evolutionary programming, and GA [18]. Each individual i maintains the following information [19]:(i)
X
i is the individual current position.
(ii)
V
i is the individual current velocity.
(iii)
Y
i is the local better position of the individual (pbest), the better position visited yet by the individual.
(iv)
Y
^ is the global better position of the swarm (gbest), the better position visited yet by the entire swarm.Figure3 shows a flowchart of the PSO algorithm.Figure 3
Flowchart of PSO algorithm.By using the above notation, the method is implemented as follows:(1)
Initialize the set constants, such as swarm size, dimension of the problem, maximum number of iterations, and upper and lower bounds.
(2)
Randomly initialize the individual positions.
(3)
Randomly initialize the individual velocities.
(4)
Repeat until the stopping condition is met.
(5)
Evaluate the fitness values using the objective function.
(6)
Determine pbest and gbest.
(7)
Determine the alteration particle velocity vector as follows:(13)
V
i
t
+
1
=
W
·
V
i
t
+
c
1
·
r
1
t
·
Y
i
t
-
X
i
t
+
c
2
·
r
2
t
·
Y
^
t
-
X
i
t
,
where,t represents current iteration, r
1
(
t
) and r
2
(
t
) represent uniform random numbers between 0 and 1, acceleration coefficients are c
1 and c
2, usually between 0 and 4, and W represents the inertia weight; a damping factor, usually decreasing from around 0.9 to around 0.4 during the computation, is calculated as follows: (14)
W
=
M
a
x
i
t
-
t
M
a
x
i
t
,
where maximum number of iterations isM
a
x
i
t.
(8)
Determine the alteration particle position vector as follows:(15)
X
i
t
+
1
=
X
i
t
+
V
i
t
+
1
.


#### 3.2.2. PSO-Based Voltage Deviation Algorithm

A flowchart of the modified PSO algorithms is shown in Figure4, and the steps are as follows.Figure 4
Flowchart of modified PSO algorithm.Step 7.
Initialize the PSO parameters, namely, swarm sizeP
s
i
z
e, dimension of the problem N
p
a
r, maximum number of iterations M
a
x
i
t, cognitive parameter c
1, social parameter c
2, upper bound value ub, lower bound value lb, and maximum velocity value V
m
a
x. Set the voltage magnitude limits, and set the possible capacitor locations, capacitor size limits, and minimum and maximum tap positions depending on the case being solved.Step 8.
Randomly create the initialized particle velocities, determine the size and location of the capacitors and tap positions of the regulators, and form particle positions as follows:(16)
p
a
r
t
i
c
l
e
i
=
C
a
p
1
,
…
,
C
a
p
m
,
l
1
,
…
,
l
m
,
T
a
p
1
,
…
,
T
a
p
n
,where m and n represent the numbers of bank capacitors and voltage regulators, respectively.Step 9.
Run a load flow using the specified load profile, perform power flow using the solution candidates, and compute the corresponding fitness value of the test system depending on the case number as in (1), (4), and (6).Step 10.
Select local best (lb) and global best (gb), and then determine alteration particle velocity vector and particle positions as in (13)–(15).Step 11.
Repeat Steps3 and 4 until the stopping condition is met.


### 3.3. Moth Flame Optimization (MFO)

#### 3.3.1. General Approach

MFO is a new population-based algorithm refined in 2015 by Mirjalili; the optimization of this algorithm reflects transverse orientation, which is the method of transmission of moths in nature at night [13]. Approximately 160,000 different groups of insects, including moths, are present in nature. Moths have two life phases: larvae and adult phases. These insects are considerably similar to the family of butterflies but possess a special feature when moving at night [20]. Moths fly straight lines over long distances by preserving a fixed angle with the moon. This mechanism is effective for traveling, especially when the light source is far. When the light source is close, moths fly around it in a spiral path and ultimately converge with it. These insects represent the candidate solutions, and their position in the search space represents the problem variables. Therefore, moths can fly in one or more dimensions by updating the position vectors. Figure 5 presents a flowchart of the MFO algorithm.Figure 5
Flowchart of MFO algorithm.This model is implemented as follows:(1)
Initialize the set constants, such as number of moths, number of variables (dimension), and upper and lower bounds.
(2)
Randomly initialize the population of moths depending on the number of moths, number of variables, and upper and lower bounds as follows:(17)
M
o
=
M
o
11
⋯
M
o
1
d
M
o
21
…
M
o
2
d
⋮
⋱
⋮
M
o
n
1
⋯
M
o
n
d
,
wheren and d represent the numbers of moths and variables, respectively.
(3)
Calculate and store the corresponding fitness values for all the moths as follows:(18)
O
M
o
=
O
M
o
1
O
M
o
2
⋮
O
M
o
n
,
wheren represents the number of moths.
(4)
Initialize the population of flames, which is equal sort population of moths, and flame fitness values, which are the equal sort moth fitness values.(19)
F
=
F
11
⋯
F
1
d
F
21
…
F
2
d
⋮
⋱
⋮
F
n
1
⋯
F
n
d
,
(20)
O
F
=
O
F
1
O
F
2
⋮
O
F
n
,
wheren and d represent the numbers of moths and variables, respectively.
(5)
Repeat until the stopping condition is met.
(6)
Calculate the distance between thejth flame and the ith moth as follows: (21)
D
i
=
F
j
-
M
o
i
.
(7)
Update the position of moths using a spiral function as follows:(22)
M
o
i
=
D
i
·
e
b
t
·
cos
⁡
2
π
t
+
F
j
,
whereD
i represents the distance, t is a random value in [
-
1,1
], and b is a constant number.
(8)
Update the flame position, which is equal to the best previous moth position and the current moth position (same as flame fitness values) as follows:(23)
F
=
S
o
r
t
M
o
i
-
1
,
M
o
i
(24)
O
F
=
S
o
r
t
O
M
o
i
-
1
,
O
M
o
i
,
wherei represents the current iteration.


#### 3.3.2. MFO-Based Voltage Deviation Algorithm

A flowchart of the modified MFO algorithms is shown in Figure6, and the steps are as follows.Figure 6
Flowchart of modified MFO algorithm.Step 12.
Initialize the MFO parameters, namely, the number of mothsN, variable number D, maximum number of repetitions M
a
x
i
t, upper bound value ub, and lower bound value lb. Set the voltage magnitude limits, and set the possible capacitor locations, capacitor size limits, and minimum and maximum tap positions depending on the case being solved.Step 13.
Randomly create the size and location of the capacitors and tap positions of the regulators, and form the initial moth position as in (17).Step 14.
Use the specified load profile to run a load flow, perform power flow using the solution candidates, and calculate the moth fitness value of the test system as in (18). Use (1), (4), and (6) depending on the case number.Step 15.
Select the best moth position as a flame position and the best moth fitness value as the flame fitness value using (23) and (24), as shown in (19) and (20), respectively.Step 16.
Calculate the distance between moths and flames, and then calculate new moth position using (21) and (22).Step 17.
Repeat Steps3–5 until the stopping condition is met.


## 3.1. Imperialist Competitive Algorithm (ICA)

### 3.1.1. General Approach

ICA was recently developed in 2007 by Esmaeil Gargari and Caro Lucas for continuous optimization problems [16]. The working philosophy corresponds to other evolutionary algorithms and initially creates random solution candidates called countries. The cost function of each solution candidate shows the power of each country. Hence, populations are composed of either colonized or imperialist countries. According to random rules, a part of a population is selected as the imperialists or the powerful countries, and the remaining part of the population comprises the colonized. Figure 1 presents a flowchart of ICA [16].Figure 1
Flowchart of ICA [16].The method is conducted as follows:(i)
Form countries: theith country is formed as follows: (8)
c
o
u
n
t
r
y
i
=
P
1
,
P
2
,
P
3
,
…
,
P
D
N
,
where DN denotes the problem variables or dimensions. Initial random values forP
i should be within the lower and upper ranges for each variable.
(ii)
Find the powers of each country by evaluating the objective function of the optimization problem as follows:(9)
f
c
o
u
n
t
r
y
i
=
f
P
1
,
P
2
,
P
3
,
…
,
P
D
N
(iii)
Select the imperialist and colonized countries. The power of a country is inversely symmetrical to its cost. The division of colonies among imperialists and the normalized value of each imperialist is defined as follows:(10)
C
n
=
c
n
-
max
⁡
c
i
,
wherec
n is the cost of nth imperialist and C
n is the normalized value.
(iv)
Then, the colony countries move to the imperialist ones to start the optimization process. The DN country population is generated, andN
i
m
p represents the most powerful population, whose members are selected as imperialists (the sets of controller coefficients with smaller cost function in this problem). The remaining N
c
o
l countries are the colonies (the sets of controller coefficients with a high cost function in this problem), each of which is a part of one of the above-mentioned empires. In the attraction policy, the colonies move toward the imperialists along Mx units and are situated in a new position. Mx is a random variable with regular distribution and can be expressed as follows: (11)
M
x
~
U
0
,
β
·
d
s
,
whereβ is a constant number greater than 1 (
β
=
2
) and ds is the space between imperialist and colony. To search different points around the imperialist, we add a random amount θ of deviation to the direction of movement as follows:
θ
~
U
-
ץ
,
ץ, where ץ is a parameter to adjust the deviation value (
ץ
=
π
/
4
).
(v)
Calculate total power of an empire. It can be determined by the power of imperialist country plus percentage of power of its colonies as follows:
T.C.n = cost (imperialist ) + Ƹ · mean (cost (colonies of empire)), where T.C.n is the total cost of nth empire and Ƹ is a positive number that is considered to be less than 1 (Ƹ
=
0.02). The weakest colony of the weakest empire is picked out.There are some other hyperparameters used in the internal calculations of this algorithm; for example, the percent of search space size is a positive number (0.02), which enables the uniting process of two empires,α is a number in the interval of 0
1 and denotes the importance of mean minimum compares to the global minimum, and revolution rate is a positive number (0.3) representing the process in which the sociopolitical characteristics of a country change suddenly.


### 3.1.2. ICA-Based Voltage Deviation Algorithm

The flowchart of the modified ICA algorithms is shown in Figure2, and the steps are as follows.Figure 2
Flowchart of modified ICA.Step 1.
Initialize the ICA parameters, namely, population sizeN
p
o
p, maximum iteration number M
a
x
i
t, number of imperialist countries N
i
m
p, and number of colony countries N
c
o
l. Set the voltage magnitude limits, and set the possible capacitor locations, capacitor size limits, and minimum and maximum tap positions depending on the case being solved.Step 2.
Randomly create the size and location of the capacitors and tap positions of the regulators, and form the initial country as follows:(12)
c
o
u
n
t
r
y
i
=
C
a
p
1
,
…
,
C
a
p
m
,
l
1
,
…
,
l
m
,
T
a
p
1
,
…
,
T
a
p
n
,where m and n represent the numbers of bank capacitors and voltage regulators, respectively.Step 3.
Run a load flow using the specified load profile and the solution candidates, perform a power flow, and calculate the fitness value of the test system depending on the case number as in (1), (4), and (6).Step 4.
Determine the imperialist and colonized countries depending on the fitness value as in (9) and (10).Step 5.
Update the size and location of the capacitors and the tap position of the regulators for all empires as in (11).Step 6.
Repeat Steps3–5 until the stopping condition is met.


## 3.1.1. General Approach

ICA was recently developed in 2007 by Esmaeil Gargari and Caro Lucas for continuous optimization problems [16]. The working philosophy corresponds to other evolutionary algorithms and initially creates random solution candidates called countries. The cost function of each solution candidate shows the power of each country. Hence, populations are composed of either colonized or imperialist countries. According to random rules, a part of a population is selected as the imperialists or the powerful countries, and the remaining part of the population comprises the colonized. Figure 1 presents a flowchart of ICA [16].Figure 1
Flowchart of ICA [16].The method is conducted as follows:(i)
Form countries: theith country is formed as follows: (8)
c
o
u
n
t
r
y
i
=
P
1
,
P
2
,
P
3
,
…
,
P
D
N
,
where DN denotes the problem variables or dimensions. Initial random values forP
i should be within the lower and upper ranges for each variable.
(ii)
Find the powers of each country by evaluating the objective function of the optimization problem as follows:(9)
f
c
o
u
n
t
r
y
i
=
f
P
1
,
P
2
,
P
3
,
…
,
P
D
N
(iii)
Select the imperialist and colonized countries. The power of a country is inversely symmetrical to its cost. The division of colonies among imperialists and the normalized value of each imperialist is defined as follows:(10)
C
n
=
c
n
-
max
⁡
c
i
,
wherec
n is the cost of nth imperialist and C
n is the normalized value.
(iv)
Then, the colony countries move to the imperialist ones to start the optimization process. The DN country population is generated, andN
i
m
p represents the most powerful population, whose members are selected as imperialists (the sets of controller coefficients with smaller cost function in this problem). The remaining N
c
o
l countries are the colonies (the sets of controller coefficients with a high cost function in this problem), each of which is a part of one of the above-mentioned empires. In the attraction policy, the colonies move toward the imperialists along Mx units and are situated in a new position. Mx is a random variable with regular distribution and can be expressed as follows: (11)
M
x
~
U
0
,
β
·
d
s
,
whereβ is a constant number greater than 1 (
β
=
2
) and ds is the space between imperialist and colony. To search different points around the imperialist, we add a random amount θ of deviation to the direction of movement as follows:
θ
~
U
-
ץ
,
ץ, where ץ is a parameter to adjust the deviation value (
ץ
=
π
/
4
).
(v)
Calculate total power of an empire. It can be determined by the power of imperialist country plus percentage of power of its colonies as follows:
T.C.n = cost (imperialist ) + Ƹ · mean (cost (colonies of empire)), where T.C.n is the total cost of nth empire and Ƹ is a positive number that is considered to be less than 1 (Ƹ
=
0.02). The weakest colony of the weakest empire is picked out.There are some other hyperparameters used in the internal calculations of this algorithm; for example, the percent of search space size is a positive number (0.02), which enables the uniting process of two empires,α is a number in the interval of 0
1 and denotes the importance of mean minimum compares to the global minimum, and revolution rate is a positive number (0.3) representing the process in which the sociopolitical characteristics of a country change suddenly.


## 3.1.2. ICA-Based Voltage Deviation Algorithm

The flowchart of the modified ICA algorithms is shown in Figure2, and the steps are as follows.Figure 2
Flowchart of modified ICA.Step 1.
Initialize the ICA parameters, namely, population sizeN
p
o
p, maximum iteration number M
a
x
i
t, number of imperialist countries N
i
m
p, and number of colony countries N
c
o
l. Set the voltage magnitude limits, and set the possible capacitor locations, capacitor size limits, and minimum and maximum tap positions depending on the case being solved.Step 2.
Randomly create the size and location of the capacitors and tap positions of the regulators, and form the initial country as follows:(12)
c
o
u
n
t
r
y
i
=
C
a
p
1
,
…
,
C
a
p
m
,
l
1
,
…
,
l
m
,
T
a
p
1
,
…
,
T
a
p
n
,where m and n represent the numbers of bank capacitors and voltage regulators, respectively.Step 3.
Run a load flow using the specified load profile and the solution candidates, perform a power flow, and calculate the fitness value of the test system depending on the case number as in (1), (4), and (6).Step 4.
Determine the imperialist and colonized countries depending on the fitness value as in (9) and (10).Step 5.
Update the size and location of the capacitors and the tap position of the regulators for all empires as in (11).Step 6.
Repeat Steps3–5 until the stopping condition is met.


## 3.2. Particle Swarm Optimization (PSO)

### 3.2.1. General Approach

PSO was originally developed in 1995 by Kennedy and Eberhart and inspired by the social behavior of schooling fish and flocking birds [17]. The birds in a group are considered an individual in the PSO method. These particles can be flown through a search space. The location of a particle in the search problem represents one solution for the problem. A new and different solution is created when the individual moves to a new location in the search space. Each solution can be evaluated using an objective function that supplies a cost of the benefit of the solution. The direction and velocity of each individual can move along all dimensions of the search space and thus can change with all generation of movement. PSO is generally considered an evolutionary computation (EC) sample. Other EC examples include evolutionary strategies, genetic programming, evolutionary programming, and GA [18]. Each individual i maintains the following information [19]:(i)
X
i is the individual current position.
(ii)
V
i is the individual current velocity.
(iii)
Y
i is the local better position of the individual (pbest), the better position visited yet by the individual.
(iv)
Y
^ is the global better position of the swarm (gbest), the better position visited yet by the entire swarm.Figure3 shows a flowchart of the PSO algorithm.Figure 3
Flowchart of PSO algorithm.By using the above notation, the method is implemented as follows:(1)
Initialize the set constants, such as swarm size, dimension of the problem, maximum number of iterations, and upper and lower bounds.
(2)
Randomly initialize the individual positions.
(3)
Randomly initialize the individual velocities.
(4)
Repeat until the stopping condition is met.
(5)
Evaluate the fitness values using the objective function.
(6)
Determine pbest and gbest.
(7)
Determine the alteration particle velocity vector as follows:(13)
V
i
t
+
1
=
W
·
V
i
t
+
c
1
·
r
1
t
·
Y
i
t
-
X
i
t
+
c
2
·
r
2
t
·
Y
^
t
-
X
i
t
,
where,t represents current iteration, r
1
(
t
) and r
2
(
t
) represent uniform random numbers between 0 and 1, acceleration coefficients are c
1 and c
2, usually between 0 and 4, and W represents the inertia weight; a damping factor, usually decreasing from around 0.9 to around 0.4 during the computation, is calculated as follows: (14)
W
=
M
a
x
i
t
-
t
M
a
x
i
t
,
where maximum number of iterations isM
a
x
i
t.
(8)
Determine the alteration particle position vector as follows:(15)
X
i
t
+
1
=
X
i
t
+
V
i
t
+
1
.


### 3.2.2. PSO-Based Voltage Deviation Algorithm

A flowchart of the modified PSO algorithms is shown in Figure4, and the steps are as follows.Figure 4
Flowchart of modified PSO algorithm.Step 7.
Initialize the PSO parameters, namely, swarm sizeP
s
i
z
e, dimension of the problem N
p
a
r, maximum number of iterations M
a
x
i
t, cognitive parameter c
1, social parameter c
2, upper bound value ub, lower bound value lb, and maximum velocity value V
m
a
x. Set the voltage magnitude limits, and set the possible capacitor locations, capacitor size limits, and minimum and maximum tap positions depending on the case being solved.Step 8.
Randomly create the initialized particle velocities, determine the size and location of the capacitors and tap positions of the regulators, and form particle positions as follows:(16)
p
a
r
t
i
c
l
e
i
=
C
a
p
1
,
…
,
C
a
p
m
,
l
1
,
…
,
l
m
,
T
a
p
1
,
…
,
T
a
p
n
,where m and n represent the numbers of bank capacitors and voltage regulators, respectively.Step 9.
Run a load flow using the specified load profile, perform power flow using the solution candidates, and compute the corresponding fitness value of the test system depending on the case number as in (1), (4), and (6).Step 10.
Select local best (lb) and global best (gb), and then determine alteration particle velocity vector and particle positions as in (13)–(15).Step 11.
Repeat Steps3 and 4 until the stopping condition is met.


## 3.2.1. General Approach

PSO was originally developed in 1995 by Kennedy and Eberhart and inspired by the social behavior of schooling fish and flocking birds [17]. The birds in a group are considered an individual in the PSO method. These particles can be flown through a search space. The location of a particle in the search problem represents one solution for the problem. A new and different solution is created when the individual moves to a new location in the search space. Each solution can be evaluated using an objective function that supplies a cost of the benefit of the solution. The direction and velocity of each individual can move along all dimensions of the search space and thus can change with all generation of movement. PSO is generally considered an evolutionary computation (EC) sample. Other EC examples include evolutionary strategies, genetic programming, evolutionary programming, and GA [18]. Each individual i maintains the following information [19]:(i)
X
i is the individual current position.
(ii)
V
i is the individual current velocity.
(iii)
Y
i is the local better position of the individual (pbest), the better position visited yet by the individual.
(iv)
Y
^ is the global better position of the swarm (gbest), the better position visited yet by the entire swarm.Figure3 shows a flowchart of the PSO algorithm.Figure 3
Flowchart of PSO algorithm.By using the above notation, the method is implemented as follows:(1)
Initialize the set constants, such as swarm size, dimension of the problem, maximum number of iterations, and upper and lower bounds.
(2)
Randomly initialize the individual positions.
(3)
Randomly initialize the individual velocities.
(4)
Repeat until the stopping condition is met.
(5)
Evaluate the fitness values using the objective function.
(6)
Determine pbest and gbest.
(7)
Determine the alteration particle velocity vector as follows:(13)
V
i
t
+
1
=
W
·
V
i
t
+
c
1
·
r
1
t
·
Y
i
t
-
X
i
t
+
c
2
·
r
2
t
·
Y
^
t
-
X
i
t
,
where,t represents current iteration, r
1
(
t
) and r
2
(
t
) represent uniform random numbers between 0 and 1, acceleration coefficients are c
1 and c
2, usually between 0 and 4, and W represents the inertia weight; a damping factor, usually decreasing from around 0.9 to around 0.4 during the computation, is calculated as follows: (14)
W
=
M
a
x
i
t
-
t
M
a
x
i
t
,
where maximum number of iterations isM
a
x
i
t.
(8)
Determine the alteration particle position vector as follows:(15)
X
i
t
+
1
=
X
i
t
+
V
i
t
+
1
.


## 3.2.2. PSO-Based Voltage Deviation Algorithm

A flowchart of the modified PSO algorithms is shown in Figure4, and the steps are as follows.Figure 4
Flowchart of modified PSO algorithm.Step 7.
Initialize the PSO parameters, namely, swarm sizeP
s
i
z
e, dimension of the problem N
p
a
r, maximum number of iterations M
a
x
i
t, cognitive parameter c
1, social parameter c
2, upper bound value ub, lower bound value lb, and maximum velocity value V
m
a
x. Set the voltage magnitude limits, and set the possible capacitor locations, capacitor size limits, and minimum and maximum tap positions depending on the case being solved.Step 8.
Randomly create the initialized particle velocities, determine the size and location of the capacitors and tap positions of the regulators, and form particle positions as follows:(16)
p
a
r
t
i
c
l
e
i
=
C
a
p
1
,
…
,
C
a
p
m
,
l
1
,
…
,
l
m
,
T
a
p
1
,
…
,
T
a
p
n
,where m and n represent the numbers of bank capacitors and voltage regulators, respectively.Step 9.
Run a load flow using the specified load profile, perform power flow using the solution candidates, and compute the corresponding fitness value of the test system depending on the case number as in (1), (4), and (6).Step 10.
Select local best (lb) and global best (gb), and then determine alteration particle velocity vector and particle positions as in (13)–(15).Step 11.
Repeat Steps3 and 4 until the stopping condition is met.


## 3.3. Moth Flame Optimization (MFO)

### 3.3.1. General Approach

MFO is a new population-based algorithm refined in 2015 by Mirjalili; the optimization of this algorithm reflects transverse orientation, which is the method of transmission of moths in nature at night [13]. Approximately 160,000 different groups of insects, including moths, are present in nature. Moths have two life phases: larvae and adult phases. These insects are considerably similar to the family of butterflies but possess a special feature when moving at night [20]. Moths fly straight lines over long distances by preserving a fixed angle with the moon. This mechanism is effective for traveling, especially when the light source is far. When the light source is close, moths fly around it in a spiral path and ultimately converge with it. These insects represent the candidate solutions, and their position in the search space represents the problem variables. Therefore, moths can fly in one or more dimensions by updating the position vectors. Figure 5 presents a flowchart of the MFO algorithm.Figure 5
Flowchart of MFO algorithm.This model is implemented as follows:(1)
Initialize the set constants, such as number of moths, number of variables (dimension), and upper and lower bounds.
(2)
Randomly initialize the population of moths depending on the number of moths, number of variables, and upper and lower bounds as follows:(17)
M
o
=
M
o
11
⋯
M
o
1
d
M
o
21
…
M
o
2
d
⋮
⋱
⋮
M
o
n
1
⋯
M
o
n
d
,
wheren and d represent the numbers of moths and variables, respectively.
(3)
Calculate and store the corresponding fitness values for all the moths as follows:(18)
O
M
o
=
O
M
o
1
O
M
o
2
⋮
O
M
o
n
,
wheren represents the number of moths.
(4)
Initialize the population of flames, which is equal sort population of moths, and flame fitness values, which are the equal sort moth fitness values.(19)
F
=
F
11
⋯
F
1
d
F
21
…
F
2
d
⋮
⋱
⋮
F
n
1
⋯
F
n
d
,
(20)
O
F
=
O
F
1
O
F
2
⋮
O
F
n
,
wheren and d represent the numbers of moths and variables, respectively.
(5)
Repeat until the stopping condition is met.
(6)
Calculate the distance between thejth flame and the ith moth as follows: (21)
D
i
=
F
j
-
M
o
i
.
(7)
Update the position of moths using a spiral function as follows:(22)
M
o
i
=
D
i
·
e
b
t
·
cos
⁡
2
π
t
+
F
j
,
whereD
i represents the distance, t is a random value in [
-
1,1
], and b is a constant number.
(8)
Update the flame position, which is equal to the best previous moth position and the current moth position (same as flame fitness values) as follows:(23)
F
=
S
o
r
t
M
o
i
-
1
,
M
o
i
(24)
O
F
=
S
o
r
t
O
M
o
i
-
1
,
O
M
o
i
,
wherei represents the current iteration.


### 3.3.2. MFO-Based Voltage Deviation Algorithm

A flowchart of the modified MFO algorithms is shown in Figure6, and the steps are as follows.Figure 6
Flowchart of modified MFO algorithm.Step 12.
Initialize the MFO parameters, namely, the number of mothsN, variable number D, maximum number of repetitions M
a
x
i
t, upper bound value ub, and lower bound value lb. Set the voltage magnitude limits, and set the possible capacitor locations, capacitor size limits, and minimum and maximum tap positions depending on the case being solved.Step 13.
Randomly create the size and location of the capacitors and tap positions of the regulators, and form the initial moth position as in (17).Step 14.
Use the specified load profile to run a load flow, perform power flow using the solution candidates, and calculate the moth fitness value of the test system as in (18). Use (1), (4), and (6) depending on the case number.Step 15.
Select the best moth position as a flame position and the best moth fitness value as the flame fitness value using (23) and (24), as shown in (19) and (20), respectively.Step 16.
Calculate the distance between moths and flames, and then calculate new moth position using (21) and (22).Step 17.
Repeat Steps3–5 until the stopping condition is met.


## 3.3.1. General Approach

MFO is a new population-based algorithm refined in 2015 by Mirjalili; the optimization of this algorithm reflects transverse orientation, which is the method of transmission of moths in nature at night [13]. Approximately 160,000 different groups of insects, including moths, are present in nature. Moths have two life phases: larvae and adult phases. These insects are considerably similar to the family of butterflies but possess a special feature when moving at night [20]. Moths fly straight lines over long distances by preserving a fixed angle with the moon. This mechanism is effective for traveling, especially when the light source is far. When the light source is close, moths fly around it in a spiral path and ultimately converge with it. These insects represent the candidate solutions, and their position in the search space represents the problem variables. Therefore, moths can fly in one or more dimensions by updating the position vectors. Figure 5 presents a flowchart of the MFO algorithm.Figure 5
Flowchart of MFO algorithm.This model is implemented as follows:(1)
Initialize the set constants, such as number of moths, number of variables (dimension), and upper and lower bounds.
(2)
Randomly initialize the population of moths depending on the number of moths, number of variables, and upper and lower bounds as follows:(17)
M
o
=
M
o
11
⋯
M
o
1
d
M
o
21
…
M
o
2
d
⋮
⋱
⋮
M
o
n
1
⋯
M
o
n
d
,
wheren and d represent the numbers of moths and variables, respectively.
(3)
Calculate and store the corresponding fitness values for all the moths as follows:(18)
O
M
o
=
O
M
o
1
O
M
o
2
⋮
O
M
o
n
,
wheren represents the number of moths.
(4)
Initialize the population of flames, which is equal sort population of moths, and flame fitness values, which are the equal sort moth fitness values.(19)
F
=
F
11
⋯
F
1
d
F
21
…
F
2
d
⋮
⋱
⋮
F
n
1
⋯
F
n
d
,
(20)
O
F
=
O
F
1
O
F
2
⋮
O
F
n
,
wheren and d represent the numbers of moths and variables, respectively.
(5)
Repeat until the stopping condition is met.
(6)
Calculate the distance between thejth flame and the ith moth as follows: (21)
D
i
=
F
j
-
M
o
i
.
(7)
Update the position of moths using a spiral function as follows:(22)
M
o
i
=
D
i
·
e
b
t
·
cos
⁡
2
π
t
+
F
j
,
whereD
i represents the distance, t is a random value in [
-
1,1
], and b is a constant number.
(8)
Update the flame position, which is equal to the best previous moth position and the current moth position (same as flame fitness values) as follows:(23)
F
=
S
o
r
t
M
o
i
-
1
,
M
o
i
(24)
O
F
=
S
o
r
t
O
M
o
i
-
1
,
O
M
o
i
,
wherei represents the current iteration.


## 3.3.2. MFO-Based Voltage Deviation Algorithm

A flowchart of the modified MFO algorithms is shown in Figure6, and the steps are as follows.Figure 6
Flowchart of modified MFO algorithm.Step 12.
Initialize the MFO parameters, namely, the number of mothsN, variable number D, maximum number of repetitions M
a
x
i
t, upper bound value ub, and lower bound value lb. Set the voltage magnitude limits, and set the possible capacitor locations, capacitor size limits, and minimum and maximum tap positions depending on the case being solved.Step 13.
Randomly create the size and location of the capacitors and tap positions of the regulators, and form the initial moth position as in (17).Step 14.
Use the specified load profile to run a load flow, perform power flow using the solution candidates, and calculate the moth fitness value of the test system as in (18). Use (1), (4), and (6) depending on the case number.Step 15.
Select the best moth position as a flame position and the best moth fitness value as the flame fitness value using (23) and (24), as shown in (19) and (20), respectively.Step 16.
Calculate the distance between moths and flames, and then calculate new moth position using (21) and (22).Step 17.
Repeat Steps3–5 until the stopping condition is met.


## 4. Experiments and Simulation Results

The proposed optimization models are experimented on IEEE 13- and IEEE 123-bus test systems. The node maps of the circuits are shown in Figures7 and 8, respectively.Figure 7
IEEE 13-bus node map.Figure 8
IEEE 123-bus node map.The different load conditions are given in Tables1 and 2 and are denoted as simulation I (minimum load) and simulation II (maximum load) on the IEEE 13- and IEEE 123-bus test systems, respectively.Table 1
Active and reactive loads on IEEE 13-bus for simulation I (minimum load) and simulation II (maximum load).
Bus number
Phases
Active load of simulation I (kW)
Active load of simulation II (kW)
Reactive load (kVar)
Load type
671
a, b, c
854
1153
660
Delta
634
a
98
160
110
Wye
634
b
79
120
90
Wye
634
c
80
120
90
Wye
645
b
106
170
125
Wye
646
b, c
160
230
132
Delta
692
a, b, c
102
170
151
Delta
675
a
320
485
190
Wye
675
b
44
68
60
Wye
675
c
202
290
212
Wye
611
c
111
169
80
Wye
652
a
80
128
86
Wye
670
a
11
17
10
Wye
670
b
42
66
38
Wye
670
c
75
117
68
WyeTable 2
Active and reactive loads on IEEE 123-bus for simulation I (minimum load) and simulation II (maximum load).
Bus number
Phases
Active load of simulation I (kw)
Active load of simulation II (kw)
Reactiveload (kVar)
Load type
Bus number
Phases
Active load of simulation I (kw)
Active load of simulation II (kw)
Reactiveload (kVar)
Load type
1
a
30
37
20
Wye
62
c
25
39
20
Wye
2
b
12
19
10
Wye
63
a
27
39
20
Wye
4
c
26
38
20
Wye
64
b
50
75
35
Wye
5
c
13
19
10
Wye
65
a
23
35
25
Delta
6
c
25
38
20
Wye
65
b
24
35
25
Delta
7
a
14
19
10
Wye
65
c
52
69
50
Delta
9
a
24
38
20
Wye
66
c
52
72
35
Wye
10
a
13
20
10
Wye
68
a
12
20
10
Wye
11
a
26
40
20
Wye
69
a
25
40
20
Wye
12
b
14
20
10
Wye
70
a
13
20
10
Wye
16
c
26
39
20
Wye
71
a
26
40
20
Wye
17
c
12
20
10
Wye
73
c
27
40
20
Wye
19
a
26
38
20
Wye
74
c
28
40
20
Wye
20
a
26
39
20
Wye
75
c
28
38
20
Wye
22
b
25
39
20
Wye
76
a
62
104
80
Delta
24
c
26
40
20
Wye
76
b
46
70
50
Delta
28
a
28
40
20
Wye
76
c
45
70
50
Delta
29
a
28
39
20
Wye
77
b
26
40
20
Wye
30
c
24
39
20
Wye
79
a
27
40
20
Wye
31
c
13
20
10
Wye
80
b
30
40
20
Wye
32
c
14
20
10
Wye
82
a
29
38
20
Wye
33
a
26
40
20
Wye
83
c
12
20
10
Wye
34
c
25
40
20
Wye
84
c
13
19
10
Wye
35
a
28
40
20
Delta
85
c
25
40
20
Wye
37
a
28
39
20
Wye
86
b
13
20
10
Wye
38
b
12
19
10
Wye
87
b
27
40
20
Wye
39
b
13
20
10
Wye
88
a
29
40
20
Wye
41
c
12
19
10
Wye
90
b
29
39
20
Wye
42
a
13
20
10
Wye
92
c
24
40
20
Wye
43
b
25
40
20
Wye
94
a
26
40
20
Wye
45
a
15
19
10
Wye
95
b
14
20
10
Wye
46
a
14
19
10
Wye
96
b
13
20
10
Wye
47
a, b, c
64
103
75
Wye
98
a
26
40
20
Wye
48
a, b, c
137
202
150
Wye
99
b
30
40
20
Wye
49
a
23
35
25
Wye
100
c
28
39
20
Wye
49
b
45
69
50
Wye
102
c
12
20
10
Wye
49
c
23
35
20
Wye
103
c
27
39
20
Wye
50
c
29
39
20
Wye
104
c
26
39
20
Wye
51
a
15
19
10
Wye
106
b
25
40
20
Wye
52
a
25
40
20
Wye
107
b
25
40
20
Wye
53
a
26
40
20
Wye
109
a
29
40
20
Wye
55
a
13
20
10
Wye
111
a
15
20
10
Wye
56
b
13
20
10
Wye
112
a
11
20
10
Wye
58
b
13
20
10
Wye
113
a
25
40
20
Wye
59
b
15
19
10
Wye
114
a
13
20
10
Wye
60
a
14
19
10
WyeGraphical representations of the bus voltage magnitudes in pu of simulations I and II with no control (test systems do not contain tap regulators and capacitor banks) are shown in Figures9 and 10 for the IEEE 13- and IEEE 123-bus test systems, respectively. The minimum and maximum voltage magnitudes in pu with no control of IEEE 13-bus test system in simulation I are 0.9081 and 0.99995, respectively, and in simulation II they are 0.89236 and 0.99993, respectively.Figure 9
Two different simulation bus voltage magnitudes of IEEE 13-bus test system with no controls.Figure 10
Two different simulation bus voltage magnitudes of IEEE 123-bus test system with no controls.The minimum and maximum voltage magnitudes in pu with no control of IEEE 123-bus test system in simulation I are 0.93317 and 0.99999, respectively, and in simulation II they are 0.91934 and 0.99999, respectively. The optimization model results for all cases, which are based on modified heuristic approaches ICA, PSO, and MFO, are graphically compared to uncontrolled results, as shown in Figures11–16 for the IEEE 13- and IEEE 123-bus test systems, respectively.Figure 11
ICA method output of IEEE 13-bus compared to uncontrolled case in simulation I condition.Figure 12
PSO method output of IEEE 13-bus compared to uncontrolled case in simulation I condition.Figure 13
MFO method output of IEEE 13-bus compared to uncontrolled case in simulation I condition.Figure 14
ICA method output of IEEE 123-bus compared to uncontrolled case in simulation II condition.Figure 15
PSO method output of IEEE 123-bus compared to uncontrolled case in simulation II condition.Figure 16
MFO method output of IEEE 123-bus compared to uncontrolled case in simulation II condition.The numerical results in Tables3 and 4 support graphically results in Figures 11–16, respectively.Table 3
Best results values of IEEE 13-bus test system in simulation I condition.
Load profile condition
Voltage magnitude in pu
Algorithm
No control
Case1 values
Case2 values
Case3 values
Simulation I (minimum load)
Minimum
No control
0.9081
0.9081
0.9081
0.9081
ICA
0.9081
0.98544
0.98462
0.98684
PSO
0.9081
0.98544
0.98477
0.97644
MFO
0.9081
0.98544
0.9846
0.9838
Maximum
No control
0.99995
0.99995
0.99995
0.99995
ICA
0.99995
1.0347
1.0209
1.0234
PSO
0.99995
1.0347
1.022
1.0275
MFO
0.99995
1.0347
1.0269
1.0227Table 4
Best results values of IEEE 123-bus test system in simulation II condition.
Load profile condition
Voltage magnitude in pu
Algorithm
No control
Case1 values
Case2 values
Case3 values
Simulation II (maximum load)
Minimum
No control
0.91934
0.91934
0.91934
0.91934
ICA
0.91934
0.97758
0.98768
0.9685
PSO
0.91934
0.97758
0.96908
0.9658
MFO
0.91934
0.97758
0.98338
0.96869
Maximum
No control
0.99999
0.99999
0.99999
0.99999
ICA
0.99999
1.035
1.0408
1.0317
PSO
0.99999
1.035
1.0404
1.0369
MFO
0.99999
1.035
1.038
1.0334The smooth curves in Figures11–16 represent the performance of Cases 1–3. Good voltage profile is observed in Case 3, in which tap regulator position and capacitor size and location are controlled. Through the curves in Figures 11–16, the voltage magnitudes can be obtained within the admissible range in any of the cases. The comparison of the simulation results that are based on ICA, PSO, and MFO is given in Figures 17–19. The numerical results in Tables 3–6 support graphically results in Figures 17–19.Figure 17
Case1 output of IEEE 13-bus compared to three methods in simulation II condition.Figure 18
Case3 output of IEEE 13-bus compared to three methods in simulation I condition.Figure 19
Case2 output of IEEE 123-bus compared to three methods in simulation I condition.The performance curves in Figures17–19 demonstrate that improved voltage is achieved in Cases 2 and 3 using ICA as shown in Tables 5 and 6. Meanwhile, as shown in Tables 5 and 6, Case 1 has the same values under all algorithms. The proposed algorithm iteration versus best fitness value in Case 3 of simulation I and simulation II of 13- and 123-bus system is shown in Figures 20 and 21, respectively. Tables 5 and 6 present comparison results, namely, best fitness values, mean voltage magnitudes for best fitness values, and standard deviation voltage magnitudes for best fitness values, in three phases under the different cases and modified heuristic approaches for the IEEE 13- and IEEE 123-bus test systems, respectively.Table 5
Best results values of IEEE 13-bus for 3 phases.
Results type
Load profile condition
Algorithm
No control
Case1 values
Case2 values
Case3 values
Best fitness values; (1) for Case 1, (4) for Case 2, (6) for Case 3
Simulation I (minimum load)
No control
0.10569
0.10569
0.10569
0.10569
ICA
0.10569
0.010363
0.006764
0.0064852
PSO
0.10569
0.010363
0.0060784
0.0068729
MFO
0.10569
0.010363
0.0056067
0.0066113
Simulation II (maximum load)
No control
0.14316
0.14316
0.14316
0.14316
ICA
0.14316
0.013653
0.00833
0.012759
PSO
0.14316
0.013653
0.0086557
0.014014
MFO
0.14316
0.013653
0.0085749
0.012991
Mean voltage magnitudes for best fitness values
Simulation I (minimum load)
No control
0.95284
0.95284
0.95284
0.95284
ICA
0.95284
1.0073
1.0101
1.0065
PSO
0.95284
1.0073
1.0063
1.0072
MFO
0.95284
1.0073
1.0048
1.0046
Simulation II (maximum load)
No control
0.94605
0.94605
0.94605
0.94605
ICA
0.94605
1.0102
1.0075
1.0067
PSO
0.94605
1.0102
1.006
1.0079
MFO
0.94605
1.0102
1.0072
1.0058
Standard deviation voltage magnitudes for best fitness values
Simulation I (minimum load)
No control
0.028612
0.028612
0.028612
0.028612
ICA
0.028612
0.015797
0.0096983
0.012119
PSO
0.028612
0.015797
0.011767
0.012222
MFO
0.028612
0.015797
0.011881
0.013154
Simulation II (maximum load)
No control
0.034842
0.034842
0.034842
0.034842
ICA
0.034842
0.017181
0.013685
0.018152
PSO
0.034842
0.017181
0.014732
0.018642
MFO
0.034842
0.017181
0.014085
0.018644Table 6
Best results values of IEEE 123-bus for 3 phases.
Results type
Load profile condition
Algorithm
No control
Case1 values
Case2 values
Case3 values
Best fitness values; (1) for Case 1, (4) for Case 2, (6) for Case 3
Simulation I (minimum load)
No control
0.39548
0.39548
0.39548
0.39548
ICA
0.39548
0.10046
0.087984
0.071935
PSO
0.39548
0.10046
0.092346
0.088997
MFO
0.39548
0.10046
0.10902
0.079843
Simulation II (maximum load)
No control
0.57663
0.57663
0.57663
0.57663
ICA
0.57663
0.1038
0.093454
0.079199
PSO
0.57663
0.1038
0.13166
0.10153
MFO
0.57663
0.1038
0.10887
0.085254
Mean voltage magnitudes for best fitness values
Simulation I (minimum load)
No control
0.96441
0.96441
0.96441
0.96441
ICA
0.96441
1.0129
1.0114
1.0089
PSO
0.96441
1.0129
1.0108
1.0113
MFO
0.96441
1.0129
1.0141
1.0124
Simulation II (maximum load)
No control
0.95697
0.95697
0.95697
0.95697
ICA
0.95697
1.0123
1.0116
1.009
PSO
0.95697
1.0123
1.016
1.0117
MFO
0.95697
1.0123
1.0122
1.0094
Standard deviation voltage magnitudes for best fitness values
Simulation I (minimum load)
No control
0.015601
0.015601
0.015601
0.015601
ICA
0.015601
0.014796
0.014398
0.014012
PSO
0.015601
0.014796
0.015417
0.014598
MFO
0.015601
0.014796
0.014815
0.012326
Simulation II (maximum load)
No control
0.018723
0.018723
0.018723
0.018723
ICA
0.018723
0.015693
0.014916
0.014876
PSO
0.018723
0.015693
0.015685
0.015844
MFO
0.018723
0.015693
0.01638
0.015424Figure 20
: The proposed algorithm iteration versus best fitness value of 13-bus test system in simulation I.Figure 21
The proposed algorithm iteration versus best fitness value of 123-bus test system in simulation II.


## 5. Conclusion

The proposed optimization model is based on three metaheuristics approaches, namely, particle swarm optimization, imperialist competitive algorithm, and moth flame optimization, for solving the voltage deviation problem. That model uses three different cases: Case1, changing the tap positions of the regulators; Case 2, changing the capacitor sizes; and Case 3, integrating Cases 1 and 2 and changing the locations of the capacitors. To prove the implementation of the proposed approach, it is applied and demonstrated on the IEEE 13- and IEEE 123-bus test systems. The numerical simulation results show that the voltage deviation problem is solved and the best solution is obtained in Case 3, which considers tap changers for the voltage regulators and the sizes and locations of the capacitors. Moreover, the ICA method provides improved results.



---

*Source: 6258350-2018-06-03.xml*