Metaheuristics and Transmission Expansion Planning: A Comparative Case Study

Abstract

Transmission expansion planning (TEP), the determination of new transmission lines to be added to an existing power network, is a key element in power system planning. Using classical optimization to define the most suitable reinforcements is the most desirable alternative. However, the extent of the under-study problems is growing, because of the uncertainties introduced by renewable generation or electric vehicles (EVs) and the larger sizes under consideration given the trends for higher renewable shares and stronger market integration. This means that classical optimization, even using efficient techniques, such as stochastic decomposition, can have issues when solving large-sized problems. This is compounded by the fact that, in many cases, it is necessary to solve a large number of instances of a problem in order to incorporate further considerations. Thus, it can be interesting to resort to metaheuristics, which can offer quick solutions at the expense of an optimality guarantee. Metaheuristics can even be combined with classical optimization to try to extract the best of both worlds. There is a vast literature that tests individual metaheuristics on specific case studies, but wide comparisons are missing. In this paper, a genetic algorithm (GA), orthogonal crossover based differential evolution (OXDE), grey wolf optimizer (GWO), moth–flame optimization (MFO), exchange market algorithm (EMA), sine cosine algorithm (SCA) optimization and imperialistic competitive algorithm (ICA) are tested and compared. The algorithms are applied to the standard test systems of IEEE 24, and 118 buses. Results indicate that, although all metaheuristics are effective, they have diverging profiles in terms of computational time and finding optimal plans for TEP.

1. Introduction

The increase in renewable power share and emerging elements, such as electric vehicles (EV) or energy storage, have created new challenges for planners and operators [1,2]. One of these key technologies is wind energy, which is growing steadily. However, wind power plants are generally built far away from demand centers, which calls for the installation of network reinforcements to be able to integrate this new production method. A transmission system is also key to optimize the use of non-dispatchable generation in general, as it can integrate systems across regions and supports efficient cross-border flows; transmission expansion planning determines the location, features, and timing of new transmission lines to supply system requirements. By considering load growth and demand forecasting, these conditions aim to minimize the investment and operational costs as well as supply electrical loads in the studied time period in a reliable and economical manner [3,4]. TEP becomes a more challenging endeavor in the presence of highly renewable shares, and, in particular, of wind energy [5,6,7], and the complexity of the problem grows [5,8].

Some recent works include, for instance, a stochastic framework incorporating wind and solar energy sources is presented in [9].

The main techniques applied to TEP are the following:

• Mathematical optimization generally finds the best expansion plan with an optimality guarantee, which metaheuristics cannot offer [10]. These methods include LP [11,12,13], NLP [14], MIP [15,16], DP [17], and Benders decomposition [18]. In this regard, Wang et al. [19] proposed a multi-stage TEP based on the deep Q-network (DQN) algorithm, by considering the construction sequence of transmission lines and using a prioritized experience replay (PER) strategy. In the suggested method, the flexibility, reliability, and economy of the expansion plans are mentioned, and the fault severity of the equipment is molded by applying the Monte Carlo method on the IEEE 24-bus reliability test system. The economic impacts of battery energy storage systems (BESS) in TEP are modeled in [20], based on using a MIP to maximize the benefits and minimize the operational restrictions of the power flow in a transmission system. The simulation case study on the Garver 6-node test system is addressed through the suggested model advantages. Han et al. [21] proposed MIP multistage stochastic programming for the TEP in the presence of uncertain demand, using dual decomposition, and by applying the branch-and-bound algorithm. The obtained results on the IEEE 30-bus network were reported and investigated. In [22], the uncertainty in different scenarios of load, fuel cost, and generation were modeled in TEP using the Kantorovich distance of social welfare distributions. The U.S. Eastern and Western Interconnections were selected as the case study, and the results were demonstrated.
• Heuristic and metaheuristic methods are based on the gradual evolution of a solution. These methods usually provide better performance and solutions of a high quality in a short computational time, albeit without an optimality guarantee. This has motivated the application of metaheuristics to complex real-world problems that cannot be solved using mathematical methods. Some specific examples that have been applied to TEP are GA [23,24], TS [25], ICA [3], GRASP [26], SA [27], swarm intelligence (SI) [28], ACO [29], DE [30], CHA [31], AIS [32], SFLA [33], SSA [34], AMOEA [35], and NSGA II [36]. Although most cases have demonstrated a good performance, in some cases, they have been shown to return local optima instead of the global optimums [37,38].

Although TEP is a design problem (and solved offline), it is important to solve it quickly because of one main reason: it is not usually solved once, but rather it is solved many times. First, as explained in [4], the usual practice involves readapting the problem many times to avoid ACLF infeasibilities. Then, the sizes of the problems to solve in the near future are expected to surge due to uncertainty (there are more scenarios) and the integration of markets (there are more nodes). As soon as a problem grows, the computation time increases exponentially. In this context, it becomes important to be time efficient; as such, there has been a wide array of literature exploiting, for instance, decomposition techniques applied to TEP problems (Reference [39] is recommended). This is why choosing an efficient method can be vital for network studies and is the motivation of this paper.

Some recent works include, for instance, a stochastic framework incorporating wind and solar energy [9]. In addition, different approaches have been proposed for modeling the wind, load, and EVs uncertainties in TEP problems. Some recent studies model the effects of wind uncertainty, as well as EVs, by applying the response program to the load [5], analyzing the effects of uncertainty in wind power and load [3,40,41,42], modeling the load correlation using unscented transformation [43], and studying the uncertainties in predicted load demand and N-K contingencies [44]. Wind and load uncertainties have been modeled using the PDF method [3,6], the successive approximation method [45], the probabilistic clustering method [37,40], and scenario reduction methods [41,44]. In addition to the effect of wind and load uncertainty, the effect of other uncertainty factors, such as the forced outage rate of transmission lines, generator, and wind speed correlation at wind farms have also been modeled [42].

Although a long list of metaheuristics has been applied to TEP in the literature, this has happened in a disparate way and often using unfair comparisons that present well-optimized methods against the most basic versions of existing ones. However, a comparison of metaheuristics in a fair setting, on exactly the same case studies, with similarly tuned parameters and on the same hardware was necessary. This comparison was missing from the literature and is the main contribution of this work. In addition, many emergent metaheuristics have not yet been tested on the TEP problem, and this paper provides a novel application for them.

In addition, a version of the problem with uncertain wind energy, demand, and EVs is considered. Then, a relatively long list of metaheuristics (GA, OXDE, GWO, MFO, EMA, SCA, and ICA) for TEP is compared. Some of these metaheuristics, such as OXDE, GWO, MFO, and EMA, have not been previously applied to TEP.

The paper is organized as follows: In Section 2, probabilistic modeling is discussed. Section 3 deals with the probabilistic objective function. Section 4 presents the applied algorithms, some of them being relatively little known. In Section 5, the results are presented. Section 6 provides an analysis of results and the final section presents the conclusions.

2. Methods/Models

This section describes the model used to model wind power, EVs, and the optimization of the TEP problem.

2.1. Power System Probabilistic Modeling

2.1.1. The Output Power of Wind Power Plant

The power generated by a wind turbine can be approximately modeled using Equation (1) [6,46]:

$$P_t=\{\begin{array}{l}0 0\le V\le V_{ci}\\ P_r\left(\frac{V-V_{ci}}{V_r-V_{ci}}\right) V_{ci}\le V\le V_r\\ P_r V_r\le V\le V_{co}\\ 0 V\ge V_{co}\end{array}$$

(1)

The wind turbine output distribution can be explained through Equation (2) [6,46]:

$$f\left(X\right)=\left\{\begin{array}{l}{F}_{zero} x=0 \\ g\left(x\right) 0\le x\le P_{rate}\\ F_{rate} x=P_{rate}\\ 0 otherwise\end{array}\right.$$

(2)

In which $g\left(x\right)$ is a fitted polynomial using discrete sampling.

2.1.2. Load Probabilistic Modeling in the Presence of EVs

The penetration of EVs at different hours of the day causes serious and adverse changes to the system load curve and leads to an extra load in the system [1]. In this case, some values are produced based on probabilistic distribution functions. Then, 10,000 samples are taken using Monte-Carlo simulation to evaluate each solution. The total average weight of the output values can be obtained using Equation (3) as the best estimate of the probabilistic load.

$$E\left(X\right)={\displaystyle \sum }_{i=1}^n{x}_i{P}_i$$

(3)

2.2. Mathematical Formulation

In the proposed method, TEP is formulated as the minimization of the sum of the investment cost in new lines (Equation (4)), which include generation costs and possibly penalties for non-supplied energy (NSE). As in most applications of TEP using optimization [4], a DCPF is applied for implementation, which has been acknowledged as a good balance between accuracy and computational complexity.

Constraints (5) to (9) add the demand balance and the DCLP laws, as well as the capacity limits for transmission and generation [3,47,48]:

$$f\left(x\right)=Min {\displaystyle \sum }_{i,j=1}^n\left(C{l}_{ij}{N}_{ij}\right)$$

(4)

$$sf+g=d$$

(5)

$$f_{ij}-{\gamma }_{ij}\left(N_{ij}^0+N_{ij}\right)\left({\theta }_i-{\theta }_j\right)=0 \forall \left(i,j\right)\in L$$

(6)

$$\left|f_{ij}\right|\le \left(N_{ij}^0+N_{ij}\right){\overline{f}}_{ij } \forall \left(i,j\right)\in L$$

(7)

$$0\le N_{ij}\le {\overline{N}}_{ij} \forall \left(i,j\right)\in L$$

(8)

$$0\le g\le \overline{g }$$

(9)

Constraint (5) presents the power balance in each node. It models Kirchhoff’s current law (KCL) in the equivalent DC network. Constraint (6) describes Ohm’s law for the equivalent DC network by considering Kirchhoff’s voltage law (KVL). All of these equations are nonlinear [3,47,48].

2.3. Proposed Algorithms

This section describes the applied algorithms, which are GA, ICA, OXDE, GWO, MFO, EMA, and SCA. GA and ICA are considered widely known, so this section will focus on the others algorithms. These descriptions will be very brief, but more information is available to interested reader in the included references.

2.3.1. OXDE

OXDE is a method based on gradual orthogonal crossover combination to improve the ability of the DE algorithm. DE was first proposed in 1995 [49,50]. This algorithm is based on the evolution of a population of solutions and includes two stages of initialization and evaluation.

The crossover operator in this algorithm is a discrete combination of a test vector $U_i\left(t\right)$ and parent vector $X_i\left(t\right)$, which produces child $X{\prime }_i\left(t\right)$, as follows:

$${X^{\prime }}_i\left(t\right)=\left\{\begin{array}{c}{U}_{ij}\left(t\right) if j\in J\\ X_{ij}\left(t\right) otherwise\end{array}\right.$$

(10)

In this equation, $X_{ij}$ refers to the i-th element of vector $X_{ij}\left(t\right)$ and J is a set of crossover points that are changed. In this method, a lookup range was first created as per Relation (11) by determining the minimum and maximum range of solutions.

$$l_{i,j}=min\left(e_i, g_i\right)+\frac{j-1}{Q-1}.\left(max\left(e_i, g_i\right)-min\left(e_i, g_i\right)\right), j=1,\dots ,Q$$

(11)

For this purpose, a combination of QOX with DE was used in the OXDE algorithm [51]. The main difference between OXDE and DE lies in its random selection, ${\overrightarrow{X}}_{i,G}$. Thus, in each generation, QOX is applied to ${\overrightarrow{X}}_{i,G}$ to produce vector $\overrightarrow{v}$. In this case, the ability to lookup for the optimal solution was strengthened.

2.3.2. GWO Algorithm

The GWO algorithm is a metaheuristic algorithm that models the hierarchical structure and social behavior of wolves at the time of hunting [52]. This algorithm is based on an initial population and four types of gray wolves $\alpha$, $\beta$, $\delta$, and $\mathsf{\omega }$, simulating the leadership hierarchy in which three steps of hunting, searching for the bait, surrounding the bait, and attacking the bait, are mimicked. Optimization is done using the movement pattern of wolves $\alpha$, $\beta$, and $\delta$. Wolf $\alpha$ is the main leader of the algorithm, while wolf $\beta$ and wolf $\delta$ are also involved and the other wolves are considered their followers. Modeling the process of surrounding the bait is done using Equations (12) and (13).

$$\overrightarrow{D}=\left|\overrightarrow{C}.{\overrightarrow{X}}_p\left(t\right)-\overrightarrow{X}\left(t\right)\right|$$

(12)

$$\overrightarrow{X}\left(t+1\right)={\overrightarrow{X}}_p\left(t\right)-\overrightarrow{A}.\overrightarrow{D}$$

(13)

where, A and C are coefficients vectors, $X_p$ is the bait location vector, $X$ is the location of each wolf, and t is the number of iterations. The vectors of coefficients A and C are calculated as follows:

$$\overrightarrow{A}=2\overrightarrow{a}.\overrightarrow{r}-\overrightarrow{a}$$

(14)

$$\overrightarrow{C}=2{\overrightarrow{r}}_2$$

(15)

Component $\mathrm{a}$ reduces linearly from 2 to 0 during successive iterations, while r and $2_{\mathrm{r}}$ are random vectors in the range of (0, 1). In fact, wolves $\alpha$, $\beta$, and $\delta$ estimate the bait location. The other wolves update their location randomly around the bait.

2.3.3. MFO Algorithm

MFO was proposed in 2015 [53] and is based on the behavior of moths around a candle light. In this algorithm, it is assumed that candidate solutions of moths are the variables of their location in a lookup space. Thus, moths can move in a multi-dimensional space by changing in their location vectors. The mechanism of updating the location of a moth around the candle is considered as a logarithmic spiral. The starting point is the moth and the ending point is the candle (as the optimal solution). This logarithmic space was considered as the following equation so that moths do not exit the lookup space while moving:

$$S\left(M_i,F_j\right)=D_i.e^{bt}.\cos \left(2\pi t\right)+F_j$$

(16)

where D shows the vector between moth I and candle j; b is a constant number to define the logarithmic spiral; and t is a random number in the range of (−1, 1). Vector D is defined as:

$$D_i=\left|F_j-M_i\right|$$

(17)

where $M_i$ shows the i-th moth and $F_j$ is the j-th candle. As value t in Equation (16) approaches –1, the moth finds a better location than the candle (the optimal solution).

2.3.4. EMA

EMA is a new optimization algorithm introduced in 2014 [54]. It models the performance of a stock exchange. The behavior of exchange market shareholders with their high, average, and low asset values is evaluated and these characteristics are used in the proposed algorithm. There are two market modes in iterations and, after each market mode, the fitness of shareholders was studied to be ranked based on their assets value. After the end of each market mode, the initial, middle and final shareholders of the population were recognized as the shareholders of group 1, group 2 and group 3. The subjects of group 1 had no tendency to exchange in any iteration while the shareholders of groups 2 and 3 exchanged stocks with separate relations. In the case of market balance mode, the algorithm absorbs shareholders to the elite, and, in the case of market fluctuations mode, the algorithm has the task of looking up for better solutions. This algorithm has two operators (the shareholders of groups 2 and 3) of the absorber and two operators of the searcher in the modes of markets including fluctuations and balance conditions. Some useful descriptions regarding this algorithm can be found in [54].

2.3.5. SCA Algorithm

The SCA algorithm was presented in 2015 [55]. In this algorithm, a mathematical model based on sin and cos functions is used for minimizing a solution. This algorithm begins by using a set of random numbers as the initial population. The location of particles is expressed by Equations (18) and (19):

$$X_i^{t+1}=X_i^t+r_1\times sin\left(r_2\right)\times \left|r_3{P}_i^t-X_i^t\right|$$

(18)

$$X_i^{t+1}=X_i^t+r_1\times cos\left(r_2\right)\times \left|r_3{P}_i^t-X_i^t\right|$$

(19)

where $X_i$ is the location of particles in the i-th dimension and j-th iteration. These two equations are combined as:

$$X_i^{t+1}=\left\{\begin{array}{c}{X}_i^{t+1}=X_i^t+r_1\times \sin \left(r_2\right)\times \left|r_3{P}_i^t-X_i^t\right|, r_4<0.5\\ X_i^{t+1}=X_i^t+r_1\times \cos \left(r_2\right)\times \left|r_3{P}_i^t-X_i^t\right|, r_4\ge 0.5\end{array}\right.$$

(20)

where $r_1$ determines the dimension of the lookup space (range of solutions); $r_2$ determines the range of movements towards the optimal solution; $r_3$ presents the weight factors of particles in moving towards the optimal point; and $r_4$ describes the amount of displacement between sin and cos in Equation (20). Since each algorithm should establish a balance between the obtained solutions for finding the optimal point, in this algorithm, Equation (21) was used to establish the balance between Equations (18) and (19).

$$r_1=a-t\frac{a}{T}$$

(21)

In which t is the i-th iteration of the algorithm, T is the total number of iterations, and $a$ is a constant.

2.4. Problem Solving Flowchart

The general flowchart of the problem-solving process for the presented algorithms is shown in Figure 1.

3. Simulation and Results

The case studies were solved using the proposed algorithms in MATLAB 7 with a computer with 3.4 GHZ processor and 1.5 GB RAM applying DCPF. It was assumed that systems will develop for the next 15 years on conditions of steadily increasing the load up to 3× the initial value in a 15-year time horizon. This very steep increase in load was assumed in order to force the problem of needing many reinforcements, which generally means a higher computational burden. This means that the results will be more interesting to extract conclusions applicable to cases with growing complexity that motivate this paper.

The system data are extracted from Reference [3]. The applied control parameter data of the proposed algorithms are presented in

Table 1. Control parameter data of the proposed algorithms.

Algorithm	System	Parameters
ICA	24-bus	zeta ($\mathsf{\zeta }$) (0.05), damp ratio (0.99), revolution rate (0.2), number of empires (8), number of countries (30), number of iterations (300)
	118-bus	zeta ($\mathsf{\zeta }$) (0.05), damp ratio (0.99), revolution rate (0.2), number of empires (5), number of countries (50), number of iterations (300)
OXDE	24 and 118-buses	population size (20), orthogonal array (Q = 3, J = 2), factor and crossover control parameter (F = CR = 0.9), number of iterations (300)
GWO	24 and 118-buses	number of search agents (30), number of iterations (300)
MFO	24-bus	number of search agents (30), number of iterations (300)
	118-bus	number of search agents (500), number of iterations (300)
EMA	24 and 118-buses	number of population (60), g1 = g2 = [0.1, 0.05], number of population in group 1 in non-oscillation mood (0.2 × number of population), number of population in group 2 in non-oscillation mood (0.3 × number of population), number of population in group 3 in non-oscillation mood (0.5 × number of population), number of population in group 1 in oscillation mood (0.2 × number of population), number of population in group 2 in oscillation mood (0.3 × number of population), number of population in group 3 in oscillation mood (0.5 × number of population)
SCA	24 and 118-buses	number of search agents (100), number of iterations (300)

. It should be noted that, in the proposed strategy, the N-1 constraint, for line and generator outages, is applied using Monte-Carlo sampling. This constraint means that no area will be islanded in the case of failures.

Tuning parameters requires accurate knowledge of the algorithm and the random search process of the problem-solving space at the time of problem-solving. In general, there is no proven and recommended method for determining the optimal values of parameters due to the random nature of the algorithms. Therefore, in this paper, in order to determine the control parameters, an iteration-based method (30 repetitions for each algorithm) and the Taguchi method [56] were applied.

3.1. The Studied Systems

In this section, the proposed algorithms were implemented on standard IEEE-RTS 24-bus and IEEE118-bus systems at the base mode without considering the load and wind uncertainty (including EVs). In these comparisons, the parameters of the proposed algorithms were selected after several iterations and trials, such as the ones shown in Table 1. The obtained results are presented in

Table 2. The results for the IEEE-RTS 24-bus system.

From	To	GA	ICA	OXDE	GWO	MFO	EMA	SCA
3	9	-	-	-	-	-	1	-
3	24	1	1	1	1	1	-	1
6	10	2	1	1	1	1	1	1
7	8	1	2	2	2	2	2	2
9	12	1	-	1	-	-	1	-
10	12	1	1	1	1	1	1	1
12	23	-	1	1	1	1	1	1
14	16	2	1	1	1	1	1	1
16	17	1	1	1	1	1	1	1
16	19		1	-	1	1	-	1
Lines		9	9	9	9	9	9	9
Cost (M$)		350	337	337	337	337	336	337
Length (Mile)		350	337	355	337	337	336	337
Time (sec)		22.3	24.9	10.40	8.90	8.47	59.28	14.79

and

Table 3. The results for the IEEE 118-bus system.

From	To	GA	ICA	OXDE	GWO	MFO	EMA	SCA
1	2	-	-	2	-	-	2	-
1	3	3	3	2	3	3	2	3
3	5	2	2	2	3	3	2	3
6	7	-	-	-	-	2	-	-
8	5	1	1	1	1	2	1	1
8	9	-	-	-	-	2	-	-
2	12	2	1	2	-	2	2	-
3	12	1	1	-	-	-	-	-
12	16	-	-	1	-	2	-	-
8	30	3	3	2	3	3	3	3
30	38	2	2	2	2	2	2	2
48	42	-	-	-	-	2	-	-
54	56	-	-	-	-	2	-	-
65	68	-	-	-	-	2	-	-
77	80	-	-	-	-	2	-	-
Lines		14	13	14	12	29	14	12
Cost (M$)		91.16	85.07	94.27	73.71	140.55	90.97	73.71
Length (Mile)		91.16	85.07	94.27	73.71	140.55	90.97	73.71
Time (sec)		491.5	473.38	55.60	52.28	794.85	206.65	225.14

. Figure 2 and Figure 3 show the convergence results of each algorithm. The total power and loads of the studied systems in a 15-year time horizon equaled 10,215 and 8550 MW for the 24-bus system and 11,148 MW and 1998.4 MW for the 118-bus system, respectively [3].

It should be noted that the investment cost of candidate transmission lines in all routes, for all cases and scenarios, is considered to be 1 M$/mile. However, in the IEEE-RTS 24-bus system, this parameter is considered to be equal to 2, 2, 4.23 and 2 M$/mile for lines 13-14, 14-23, 16-23 and 19-23, respectively [3].

In Table 2, regarding the optimal results for the 24-bus network, all of the algorithms resulted in the same optimal solution according to the network conditions (this case is provided as a criterion for the initial tests of the proposed algorithms). According to these results, the initial parameters of all algorithms were well adjusted and the same optimal results were proposed in terms of the number of candidate lines, their locations and investments costs. Only in the obtained results from the OXDE algorithm was there a difference with other algorithms in terms of the location of one of the new candidate lines, but it was the same in terms of investment costs. These results showed the ability of all the proposed algorithms to solve the TEP problem in the initial network conditions. However, this trend is different in Table 3, for the 118-bus network, due to the larger network size and the increase in the number of decision variables of the algorithms.

As shown in Table 2 and Table 3, the results indicate that all proposed algorithms had a good performance in solving the TEP discrete problem in small-sized networks. In comparison with the ICA algorithm in Reference [3], where the results were compared with the results of other metaheuristic algorithms, similar results were achieved. However, the GWO and SCA algorithms, and at the next stage OXDE algorithm, had better performances when applied to larger and practical networks; the reason for this was the better global search of these algorithms.

3.2. Considered Scenarios

The cases were studied for the following three scenarios.

• First scenario: The corrected system considering uncertainty in wind power plants.
• Second scenario: The corrected system considering uncertainty in the loads connected to buses without generators due to load uncertainty, and EVs equal to 10% of the total loads connected to the system.
• Third scenario: The corrected system considering uncertainty in the wind energy, load, and EVs (in non-generation buses).

Then, the cases were implemented on the 24- and 118-bus systems to study the effects of different operating conditions in addition to the system’s base mode. The corrected IEEE-RTS 24-bus system is based on Figure 4. This system was corrected in buses 1 and 15, adding 2 wind farms with capacities of 180 MW (with characteristics of $V_{ci}=4$, $V_r=10$, and $V_{co}=22$ m/s at average wind speed of 5.4 m/s) [6]. The corrected IEEE 118-bus system was changed for buses 36, 69, and 77, adding three wind farms with capacities of 250, 180, and 100 MW, respectively [3,57]. In these systems it was assumed that EVs in non-generation buses are connected and that the charges equaled 10% of the total system loads.

3.3. Simulation and Results

The results obtained from optimizing the proposed algorithms are categorized and presented in the studied scenarios on the 24- and 118-bus systems as follows:

3.3.1. First Scenario

The results obtained from this scenario, applying the proposed algorithms, are presented in

Table 4. The results for the IEEE RTS 24-bus system in the first scenario.

From	To	GA	ICA	OXDE	GWO	MFO	EMA	SCA
1	5	1	1	1	1	1	1	1
3	24	1	1	1	1	1	1	1
6	10	1	1	1	1	1	1	1
7	8	2	2	2	2	2	2	2
10	12	1	1	1	1	1	1	1
14	16	1	1	1	1	1	1	1
15	24	1	1	1	1	1	1	1
16	17	1	1	1	1	1	1	1
Lines		9	9	9	9	9	9	9
Cost (M$)		332	332	332	332	332	332	332
Length (Mile)		332	332	332	332	332	332	332
Time (sec)		20.56	25.99	8.54	9.08	10.54	74.72	35.58

and

Table 5. The results for the IEEE 118-bus system in the first scenario.

From	To	GA	ICA	OXDE	GWO	MFO	EMA	SCA
1	2	-	-	1	-	1	1	-
1	3	3	3	2	3	2	2	3
3	5	2	2	2	2	2	2	3
4	5	-	-	-	-	2	-	-
8	5	1	1	1	1	1	1	1
9	10	-	-	-	-	2	-	-
2	12	2	1	1	1	2	1	-
3	12	1	1	-	-	-	-	-
23	24	-	-	-	-	-	1	-
8	30	3	3	2	2	2	2	2
30	38	3	2	2	2	2	1	2
34	36	-	-	-	-	2	-	-
34	37	-	-	-	-	2	-	-
60	61	-	-	-	-	2	-	-
64	61	-	-	-	-	2	-	-
77	78	-	-	-	-	2	-	-
108	109	-	-	-	-	2	-	-
114	115	-	-	-	-	2	-	-
Lines		15	13	11	11	31	11	11
Cost (M$)		96.56	85.07	69.78	64.03	112.476	69.03	68.67
Length (Mile)		96.56	85.07	69.78	64.03	112.476	69.03	68.67
Time (sec)		325.9	662.02	93.88	50.95	863.93	175.79	223.55

. In addition, Figure 5 and Figure 6 depict the convergence results of each algorithm to solve the first scenario.

3.3.2. Second Scenario

The results obtained from this scenario, applying the proposed algorithms, are presented in

Table 6. The results for the IEEE-RTS 24-bus system in the second scenario.

From	To	GA	ICA	OXDE	GWO	MFO	EMA	SCA
1	2	1	1	1	1	1	1	1
6	10	1	1	1	1	1	1	1
7	8	3	2	2	2	2	2	2
10	12	1	1	1	1	1	1	1
12	13	1	1	1	1	1	1	1
14	16	1	1	1	1	1	1	1
16	17	1	1	1	1	1	1	1
Lines		9	8	8	8	8	8	8
Cost (M$)		273	257	257	257	257	257	257
Length (Mile)		273	257	257	257	257	257	257
Time (sec)		24.32	21.65	8.31	15.53	10.80	47.38	23.44

and

Table 7. The results for the IEEE 118-bus system in the second scenario.

From	To	GA	ICA	OXDE	GWO	MFO	EMA	SCA
1	2	-	-	2	-	1	2	-
1	3	3	3	2	3	2	2	3
3	5	2	2	2	2	2	2	3
8	5	1	1	1	1	1	1	1
9	10	-	-	-	-	2	-	-
11	12	-	-	-	1	2	-	-
2	12	1	-	2	1	3	1	-
3	12	1	1	-	-	-	-	-
16	17	1	1	-	-	-	-	-
8	30	3	2	2	2	2	2	3
30	38	2	2	2	2	2	2	2
29	31	-	-	-	-	2	-	-
34	36	-	-	1	-	-	-	-
35	36	-	-	1	-	-	-	-
34	37	-	-	-	-	2	-	-
48	49	-	-	-	-	1	-	-
54	56	-	-	-	-	2	-	-
63	59	-	-	-	-	2	-	-
63	64	-	-	-	-	2	-	-
68	69	-	-	-	-	2	-	-
108	109	-	-	-	-	1	-	-
114	115	-	-	-	-	2	-	-
68	116	-	-	-	-	2	-	-
Lines		14	12	15	12	35	12	12
Cost (M$)		97.56	91.88	75.94	65.99	132.81	79.77	73.71
Length (Mile)		97.56	91.88	89.63	65.99	132.81	79.77	73.71
Time (sec)		311.4	933.74	149.10	105.74	721.44	173.25	184.44

. In addition, Figure 7 and Figure 8 depict the convergence results of each algorithm to solve the second scenario.

3.3.3. Third Scenario

The results obtained from this scenario, applying the proposed algorithms, are presented in

Table 8. The results for the IEEE-RTS 24-bus system in the third scenario.

From	To	GA	ICA	OXDE	GWO	MFO	EMA	SCA
1	5	1	1	1	1	1	1	1
6	10	1	1	1	1	1	1	1
7	8	2	2	2	2	2	2	2
9	11	1	1	1	1	1	1	1
10	12	1	1	1	1	1	1	1
14	16	1	1	1	1	1	1	1
16	17	1	1	1	1	1	1	1
Lines		8	8	8	8	8	8	8
Cost (M$)		260	260	260	260	260	260	260
Length (Mile)		260	260	260	260	260	260	260
Time (sec)		24.4	24.4	8.44	9.80	8.38	57.08	21.50

and

Table 9. The results for the IEEE 118-bus system in the third scenario.

From	To	GA	ICA	OXDE	GWO	MFO	EMA	SCA
1	2	-	-	1	-	1	1	-
1	3	3	3	2	3	2	2	3
3	5	2	2	1	2	1	1	2
4	5	-	-	-	-	2	-	-
8	5	1	1	-	1	-	-	1
8	9	-	-	-	-	1	-	-
9	10	-	-	-	-	2	-	-
2	12	-	-	1	-	2	1	-
3	12	1	1	1	-	-	1	-
12	14	2	-	2	-	-	1	-
12	16	1	-	1	-	-	-	-
14	15	-	-	-	-	2	-	-
15	17	-	-	-	-	1	-	-
16	17	-	1	-	-		-	-
8	30	1	2	2	2	2	2	2
30	38	2	2	1	1	1	1	1
35	36	-	-	-	-	3	-	-
54	56	-	-	-	-	2	-	-
55	56	-	-	-	-	1	-	-
60	61	-	-	-	-	1	-	-
63	64	-	-	-	-	2	-	-
68	69	-	-	-	-	2	-	-
68	81	-	-	-	-	1	-	-
108	109	-	-	-	-	2	-	-
114	115	-	-	-	-	2	-	-
68	116	-	-	-	-	2	-	-
Lines		13	12	12	9	35	10	9
Cost (M$)		93.43	91.88	73.98	52.47	141.426	73.98	52.47
Length (Mile)		93.43	91.88	89.39	52.47	141.426	73.98	52.47
Time (sec)		389.2	656.89	96.19	49.99	849.88	173.85	205.30

. In addition, Figure 9 and Figure 10 depict the convergence results of each algorithm to solve the third scenario.

4. Empirical Results

The performances of the proposed algorithms are categorized as follows:

4.1. The Number of Iterations

As metaheuristics are affected by the random evolution of an algorithm, case studies have been solved for a series of simulations. For this purpose, each optimization algorithm was implemented, separately, several times. Based on the results, it can be stated that OXDE, GWO, and EMA obtained the optimal solution in at least five separate iterations. However, the MFO algorithm, and in particular SCA, obtained the optimal solution in at least 10–15 separate iterations.

4.2. Convergence in Finding Optimal Solution

In the IEEE 24-bus system, it can be stated that OXDE, GWO, MFO, and EMA algorithms converged to the optimal solution. However, in the 118-bus system, this was not true. Only GWO and SCA converged to the optimal solution. This is depicted in Figure 11 and Figure 12.

4.3. In Terms of Implementation Time

Figure 13 shows the implementation time related to the 24-bus system and Figure 14 shows the 118-bus system. The GWO and OXDE algorithms have higher speeds in obtaining optimal solutions (while being applied to different systems) than the other proposed algorithms.

4.4. Wilcoxon Signed-Rank Test

Some statistical analyses were performed on the obtained results using simulations to measure the performance of different algorithms, OXDE, GWO, MFO, EMA and SCA. For this purpose, the Wilcoxon signed-rank test was used. This test is used to examine two dependent samples or to match two samples to determine the significant differences of different algorithms. If an algorithm establishes a statistical significance by providing sufficient evidence against the null hypothesis, then it will be considered a strong algorithm [58]. In the Wilcoxon signed-rank test, a probability (p-value) less than 0.05 indicates the reality of the observed difference and is considered as a difference due to chance, randomness, or sample error. The statistical significance shows the extent to which statistical findings may be the result of chance and accident. If the data are in the range of 95% of the samples and the null hypothesis is correct, then the statistical samples are not significant at the 5% level and the null hypothesis is not rejected. If the data sample are at the outer 5% range of the samples, then the sample does not seem to support the null hypothesis. In this case, the opposite hypothesis is usually accepted, the null hypothesis is rejected, and the statistical findings are meaningful [58,59].

Figure 15 depicts this concept, in which the blue line is the probability distribution for the test statistics, the critical region is the set of values to the right of the observed data point and the p-value is highlighted.

Table 10. The results of the p-value in the Wilcoxon signed-rank test.

Test System	Scenario	ICA	OXDE	GWO	MFO	EMA	SCA
IEEE 24-bus system	1	1.272 × 10−64	4.7131 × 10−55	2.4524 × 10−58	7.8016 × 10−54	3.3153 × 10−59	6.7927 × 10−54
	2	4.1495 × 10−52	4.3037 × 10−56	1.2041 × 10−58	5.1854 × 10−51	1.6555 × 10−52	2.1612 × 10−51
	3	4.1875 × 10−52	1.118 × 10−56	6.4352 × 10−58	2.266 × 10−51	9.1537 × 10−56	1.9928 × 10−51
IEEE 118-bus system	1	7.0049 × 10−52	4.1162 × 10−51	5.4061 × 10−53	6.0229 × 10−51	5.2716 × 10−52	2.0449 × 10−54
	2	8.8571 × 10−56	5.1188 × 10−51	2.934 × 10−52	6.0497 × 10−51	3.8101 × 10−51	8.3324 × 10−58
	3	1.5863 × 10−55	2.4559 × 10−51	9.7787 × 10−54	6.0266 × 10−51	2.7647 × 10−52	8.6361 × 10−53

presents the p-values obtained using the Wilcoxon signed-rank test with a marked rank. As is clear from the results, the p-value obtained was much less than the desired value of 0.05 and proves the statistical significance of the obtained results.

5. Discussion

This paper compared the performance of several metaheuristics, OXDE, GWO, MFO, EMA, SCA, and ICA [3], in solving TEP problem in the presence of uncertainty in power generation by wind power plants, as well as demand and EVs, as one of the consumers with higher uncertainty.

All the metaheuristics studied were able to find the optimal solution as identified by classical optimization (MIP). According to the obtained results, the performances of GWA and OXDE were remarkably robust. The exploration factor indicated the generality of the search area and the extraction parameter described the ability to find better solutions around a specific solution. The exploration operation searches for new, better solutions. Extraction follows small and important changes to the current solution. This theorem was well defended in the case of the OXDE algorithm for the two case studies, using the IEEE 24 and 118 buses, in which the algorithm was able to achieve higher-quality solutions, regardless of the number of problem variables. The remaining algorithms were considerably slower across the case studies, with ICA and EMA performing considerably worse. Of course, the important factor about metaheuristic algorithms is that they do not provide completely definite results for use, an issue that needs to be improved by other researchers.

6. Conclusions

The exercise carried out in this paper has allowed to compare several metaheuristic algorithms when solving TEP in terms of variability, convergence, and computation time. Most metaheuristics were acceptable in terms of both total cost and computation time, with some exceptions, such as ICA and MFO, which performed worse in the large case studies. The rest of the algorithms performed acceptably well, with OXDE, for instance, displaying a robust performance across case studies, consistently avoiding local optima. Metaheuristics can be a good additional tool in the arsenal for solving TEP problems, even when combined with mixed-integer programming or heuristic methods. When solving complex problems, having a diversity of tools is key, and considering metaheuristics increases the options for solving the TEP problems effectively and efficiently.