Hybrid Imperialist Competitive and Grey Wolf Algorithm to Solve Multiobjective Optimal Power Flow with Wind and Solar Units

The optimal power flow (OPF) module optimizes the generation, transmission, and distribution of electric power without disrupting network power flow, operating limits, or constraints. Similarly to any power flow analysis technique, OPF also allows the determination of system’s state of operation, that is, the injected power, current, and voltage throughout the electric power system. In this context, there is a large range of OPF problems and different approaches to solve them. Moreover, the nature of OPF is evolving due to renewable energy integration and recent flexibility in power grids. This paper presents an original hybrid imperialist competitive and grey wolf algorithm (HIC-GWA) to solve twelve different study cases of simple and multiobjective OPF problems for modern power systems, including wind and photovoltaic power generators. The performance capabilities and potential of the proposed metaheuristic are presented, illustrating the applicability of the approach, and analyzed on two test systems: the IEEE 30 bus and IEEE 118 bus power systems. Sensitivity analysis has been performed on this approach to prove the robustness of the method. Obtained results are analyzed and compared with recently published OPF solutions. The proposed metaheuristic is more efficient and provides much better optimal solutions.


Introduction
Optimal power flow (OPF) is the mathematical tool used to find the optimal settings of the power system network [1]. The main target of the OPF problem is to optimize a specific objective function while satisfying feasibility and security constraints [2]. OPF has been broadly used in previous studies [3], and has served as a substantial optimization test problem because it is characterized as multidimensional, large-scale nonlinear nonconvex, and highly constrained [4,5]. Several OPF formulations have been developed during the last few decades in order to optimize the operation of an electric power system subject to physical constraints [6]. The emerging optimization problem uses different names and different objective functions [7]. A lot of OPF solution approaches have been developed, each with distinct mathematical characteristics and computational requirements [8,9]. In recent years, OPF optimization problems have regained importance due to the rapid adoption of distributed energy resources in the network [10]. The integration of distributed and intermittent renewable energy sources, such as photovoltaic (PV) systems and wind energy (WE), into modern power systems has introduced new types of challenges for operating and managing the power grid [11]. The stochastic nature of WE and PV units must be taken into consideration to ensure successful implementation of these intermittent energy sources to the network [12]. Solving the OPF problem has become more complicated with massive incorporation of renewable resources that impose volatile dynamics to the power grid because of their uncertainty.
Conventional optimization methods, like linear (LP) and nonlinear programming (NLP) [13], quadratic programming (QP) [14], interior point method (IPM), and Newton's method [15] showed excellent convergence characteristics in solving OPF problems; however, they use theoretical assumptions not suitable for practical systems having non-differentiable, non-smooth, and nonconvex objective functions. Sometimes, the preceding approaches fail to represent the main characteristics of the fuel cost as a convex function [16]. Such a situation emerges when piecewise quadratic cost, valve points, and prohibited operating zones characteristics are presented [17]. Usually, multiple trials and accurate adjustment of associated parameters are needed to achieve the optimal solution for a specific problem. As a result, we need a faster and more robust algorithm to solve realistic OPF problems. Recently, many publications have focused on metaheuristics to solve hard optimization problems. Metaheuristics, based on a common set of principles which make it possible to design solution algorithms, may be used to overcome the abovementioned weaknesses. Most metaheuristics have the following features: they are inspired from nature, they do not use the objective function's Hessian or gradient matrix, they make use of stochastic components, and they have many parameters that need to be adapted to the problem [18]. The following artificial intelligence based optimization methods have been successfully used to solve OPF problems: moth swarm algorithm, MSA [19]; modified particle swarm optimization, MPSO [20]; modified differential evolution, MDE [21]; moth-flame optimization, MFO [22]; flower pollination algorithm, FPA [23]; adaptive real coded biogeography-based optimization, ARCBO and real coded biogeography-based optimization, RCBBO [24]; grey wolf algorithm, GWO and differential evolution, DE [25]; modified Gaussian bare bones imperialist competitive algorithm, MGBICA and Gaussian bare bones imperialist competitive algorithm, GBICA [26]; artificial bee colony, ABC [27]; simulated annealing and hybrid shuffle frog leaping algorithm [28]; Lévy mutation teaching-learning-based optimization, LTLBO [29]; teaching learning-based optimization, TLBO [30]; hybrid MPSO and shuffle frog leaping algorithms, HMPSOSFLA, and particle swarm optimization, PSO [31]; Gbest-guided artificial bee colony, GABC [32]; differential search algorithm, DSA [33]; efficient evolutionary algorithm, EEA and eclectic genetic algorithm, EGA [34]; particle swarm optimization with aging leader and challengers, ALCPSO [35]. The above optimization approaches have been developed to solve simple and multiobjective OPF problems. These algorithms performed better than traditional mathematical programming techniques in solving multiobjective optimization problems because they are less affected by the Pareto front shape, and are capable of finding the optimal solutions sets in one run [36]. The assessment of these metaheuristics is commonly based on experimental comparisons.
The objective of this research is to develop an original metaheuristic called hybrid imperialist competitive and grey wolf algorithm (HIC-GWA) to solve twelve different cases of simple and multiobjective OPF problems for hybrid power systems that includes PV and WE sources, in order to find effective, faster, and better solutions. The potential and efficiency of the HIC-GWA are presented and evaluated on two standard test systems: IEEE 30 and IEEE 118 bus power systems. Simulation results are compared with the abovementioned optimization approaches. The proposed HIC-GWA is a combination of two algorithms: the imperialist competitive algorithm (ICA) and the grey wolf optimization (GWO). ICA is a sociopolitically inspired optimization strategy that has been proposed to handle tough optimization problems [37]. This approach exhibits good performance in terms of convergence rate and improved global optimum [38,39]. The GWO algorithm is an original swarm intelligence technique stimulated by the leadership hierarchy and hunting structure of grey wolves. This robust algorithm has been used in different complex problems because of the reduced number of random parameters and a faster convergence due to continuous reduction of search space [40,41]. Each optimization technique, ICA and GWO, possesses certain specific intelligence to search for the solution of a problem. Therefore, a collection of such abilities enhances the power of the proposed metaheuristic.

Objective Functions
OPF research seeks to compute a steady state operating point that reduces cost, emission, loss, etc., while maintaining good system performance. The general OPF problem usually contains discrete and continuous control variables. It is a large-scale, nonconvex, and nonlinear optimization problem. OPF seeks to optimize the generation, transmission, and distribution of electric power with no disruption of flow, operating limits, or constraints. Similar, to other power flow analysis techniques, OPF also allows the determination of system's state of operation, that is, the injected power, voltage, and current throughout the electric power system. In this context, a large array of OPF formulations and solution methods have been developed. Furthermore, OPF research is growing, due to contemporary electricity markets and integration of renewable energy sources.
The following objective functions are minimized by the proposed HIC-GWA: Wind energy is increasingly being integrated into the power grid due to its rapidly declining cost and emission free nature. The WE power cost function can be modeled as Wind power operators get penalized if they fail to provide the scheduled amount of wind energy. Penalty costs consists of two parts: (1) underestimation cost which should be considered when available power of wind farm is not utilized, (2) overestimation cost which is calculated for buying power from alternate sources (reserves) or load shedding. These costs can be modeled as follows [12]: where i = 1, 2, . . . , n w and f (P) symbolizes the probability density function (PDF) of WE output power. The WE total cost is given by the following function: To model the unpredictable nature of wind speed, we use the Weibull distribution with PDF f (V w ) and cumulative distribution function (CDF), F(V w ), defined as follows [12]: The generated power of WE is computed as where Energies 2018, 11, 2891 4 of 23 V w and V w,r symbolizes speed and rated speed of WE generators, V w,in and V w,out symbolizes cut-in and cut-out speed of WE generators, K, C symbolizes shape and scale parameters of the Weibull distribution.

PV Cost Function
Photovoltaic systems are gaining popularity as a clean energy source due to their affordable cost and simple design. PV characteristics are highly dependent on various factors, including irradiance level, shades, and temperature, which makes it hard to accurately forecast its power production. The generation and penalty costs for PV power can be calculated as follows: C ue,pv,i = K ue,pv,i P pv,r,i P pv,i P − P pv,i f (P)dP (9) C oe,pv,i = K oe,pv,i P pv,i 0 P pv,i − P f (P)dP (10) where i = 1, . . . , n v and f (P) represent the PDF of the PV unit's output power.
The total cost of PVs is given by the following function: The PDF of the ith PVs' output power is calculated as follows: • Solar cells or PV cells are hypersensitive to the amount of solar radiation. The PDF of solar radiation f (R) can be modeled by a beta distribution [12]: where Γ(.) is the gamma function, α and β are parameters of the beta distribution, and R is the solar radiation.

•
The relation between power output of PV and output power of solar cell generator which is related to the solar radiation can be calculated as follows: where R C and R STD are solar radiation in W/m 2 . Usually, a typical solar radiation point is set to 150 W/m 2 , and it is set to 100 W/m 2 under standard conditions.

Basic Fuel Cost Function
The basic fuel cost is OPF's most common objective function. A power plant's fuel cost is commonly modeled as a quadratic function [42]: where i represents the ith power plant and n G is the number of power plants. a i , b i , and c i are cost coefficients for the ith power plant. P Gi is power of ith power plant.

. Piecewise Quadratic Fuel Cost Function
For a given operating range, power plants usually use the most economical available fuel option. Such a system has piecewise quadratic fuel cost function Each quadratic piece of the fossil fuel cost can be calculated using the following function: where n f is the number of fossil fuel options for ith power plant and a i,k , b i,k , c i,k , are coefficients for the cost of ith power plant for kth fuel option.

Piecewise Quadratic Fuel Cost with Valve Point Loading
The generator cost is a convex function with an incremental heat rate curve, subjected to discontinuities caused by the steam admission valves in large turbines. The valve point effect must be included in order to have an accurate cost for each generating unit [43]: where e i and f i are valve point cost coefficients of ith power plant.

Emission Cost Function
To produce electricity, a fossil fuel power station burns natural gas, petroleum, or coal. Significant amounts of emission are produced during the burning process. In this paper, the emission level of the two important pollutants, nitrogen oxides (NOx) and sulfur oxides (SOx), are modeled by the following function [19]: where, α i , β i , ζ i , and θ i are emission coefficients of ith power plant.

Power Loss Cost Function
To reduce the active power loss of transmission lines, the following power loss function has to be minimized [27]: where n l is the number of transmission lines, (G ij ,B ij ) are (real, imaginary) of ith jth components of the admittance matrix, δ ij is the angle separating the ith bus from the jth bus, and V i is the ith bus voltage.

Fuel Cost and Active Power Loss Cost Function
This function model two simple objectives: fuel cost and active power loss.
where β 1 is a weighting factor. One of the valuable quality and security indices is the voltage magnitude fluctuation from the specified reference value at each load bus. This function models both fuel cost and voltage deviation (VD).
where n L is the number of load buses, V Li is the ith voltage of load buses, and β 2 is a weighting factor.

Fuel Cost and Voltage Stability Enhancement
Voltage stability is the ability of a power system to sustain stable voltages at each bus within acceptable level after being exposed to a disruption. It is represented by indices like the L index, which has been introduced to evaluate the stability limit [19]. The L index is a quantitative measure of how close a point is to the system stability limit. Reducing the value of the L index is very important in power system planning and operations.
This function models the fuel cost and the L index maximum.
where β 3 is a weighting factor. The nodal admittance relates system voltages and currents as Equation (23) can be reformulated by separating the PQ bus-active and reactive power; and the PV bus-active power and voltage magnitude.
The L index is calculated by where Y_1 and Y_2 are the system Y bus submatrices.

Fuel Cost and Voltage Stability Enhancement during Contingency Condition
Transmission lines outages are used to replicate a contingency condition. This function models both fuel cost and enhancement of voltage stability.
where β 4 is a weighting factor. This function models fuel cost, emission, voltage deviation, and active power loss.

Constraints
The OPF optimization problem should satisfy the following constraints: (1) Active and reactive power balances where the number of power system bus is represented by n. P Gi , Q Gi , and P Di , Q Di are active and reactive power of generators and load, respectively, at the ith bus.
(2) The voltage magnitude of the power plant where V min i and V max i are minimum and maximum limit of ith bus voltage of power plants V i .

(3) Prohibited operating zones
There is a risk of machine or accessory failure when a power plant operates outside acceptable ranges, as shown in Equations (32)- (41).
where P l Gi,k and P u Gi,k are lower and upper bounds of the kth POZ of ith unit. P min Gi and P max Gi are active power boundaries of ith generator.
(4) Active and reactive power where Q min Gi and Q max Gi are boundaries' reactive power of ith traditional generator.
where Q min c,i and Q max c,i are the ith shunt compensator Q c,i limits. N cap represents the number of capacitors linked to the network.
where S max i is MVA's maximum. N l is the number of lines.
Each wind turbine is equipped with a squirrel cage induction generator modeled as PQ buses [44].
where X i is the sum of the leakage reactance of the stator and rotor of the ith wind turbine. V ww,i and Q w,i represents the voltage and the reactive power of the associated bus of the ith wind generator.

New Hybrid Optimization Algorithm
In this research, a new metaheuristic HIC-GWA is considered to solve twelve cases of simple and multiobjective OPF problems. This approach is a combination of two algorithms: ICA and GWO. Each of such optimization techniques, ICA and GWO, possesses certain specific heuristics to search for the solution of a problem. Therefore, a collection of such abilities enhances the power of the proposed metaheuristic.

Imperialist Competitive Algorithm (ICA)
The ICA is stimulated by the sociopolitical aspect of imperialistic competition between countries in the same population. Countries can be colonies or imperialists. Powerful countries are selected to be imperialists. Colonies are distributed among imperialists based on imperialist's power. Empires are formed with imperialist states and their colonies. Imperialistic competition between empires converge to one imperialist state which represent the optimum point of the ICA [37][38][39].

Creation of Initial Empires
A country is usually represented by an N var -dimensional array of variables that should be optimized. country = [P 1 , P 2 , . . . , P N var ] The cost of each country is inversely proportional to its power.
The cost function f is given by In the initialization process, the algorithm produces N Country initial countries. A certain number of empires, N imp , are formed with the most powerful countries. The remaining countries, N col , become colonies of the empires.
The cost of the nth imperialist is The power of the nth imperialist is The nth empire's initial number of colonies is where N col is the total number of original colonies.

Assimilation
To absorb their colonies, the imperialist states use different sociopolitical axes to make colonies move toward themselves. This movement can be modeled using different optimization axes. In a two-dimensional problem, colonies are absorbed by the imperialist using language and culture. Colonies will move toward the imperialist among these two axes. This acclimatization, modeled by approaching the colonies to the imperialist, will continue until all colonies are fully assimilated. This motion is represented by a uniform distribution: where β > 1. d represents the distance separating the colony to the imperialist state. A random deviation θ is added to the direction of movement to increase the search space around the imperialist. θ is represented by a uniform distribution.
where γ accommodates the fluctuation from the initial direction.

Revolution
Revolution is simulated to denote a shift in sociopolitical institutions that prohibits the convergence of a country to a local minimum which increases the exploration of this approach.

Exchanging Positions of a Colony and the Imperialist
The colony and the imperialist countries will change positions if the colony reaches a position with higher power than the imperialist.

Union of Empires
While moving toward the optimum solution, two imperialists may merge into one empire if they are too close to each other. Their colonies become colonies of the new empire which take the position of one of the two imperialists.

Total Empire Power
An empire's total power is highly correlated to the power of the imperialist country, but it is slightly affected by the power of the colonies. An empire's total cost is modeled as where ξ is a positive small weight factor.

Imperialistic Competition
This competition is built on the total power of the empires. Empires try to take control of each other's colonies to expand their territory. Every empire will have the possibility of possessing colonies that it is competing for. Powerful empires will control weaker colonies. The weakest colony of the weakest empire will be selected in the initiation process of the competition. An empire's possession probability (PP) is proportional to the empire's total power.
Empire's normal total cost: Empire's possession probability: The algorithm will stop after a predetermined number of iterations which represents maximum number of decades.

Grey Wolf Optimizer (GWO)
The GWO is a conventional swarm intelligence algorithm stimulated by the leadership hierarchy and hunting structure of grey wolves. This algorithm has been used in diverse complex problems because of its simplicity and robustness. The wolf colony (N w ) is divided into four clusters: alpha (α), beta (β), delta (δ), and omega (Ω). The hunting mechanism involves three main steps: searching and approaching the prey, encircling and harassing the prey, and attacking the prey [40,41].

Social Hierarchy
The leaders α are mostly responsible for making decisions about hunting. They are considered as the fittest solution. The second-best candidates are the β wolves, based on the democratic behavior of the colony. Consequently, the δ wolves take place after the β wolves. The rest are assumed to be the ωwolves. The optimization (hunting) process is guided by α, β, and δ, with the ω wolves tracking them.

Encircling Prey
Hunting in groups is another compelling social behavior of grey wolves. A grey wolf can revise its position neighboring the prey in any random place using the following equations [40]: where → r 1 and → r 2 are random vectors in [0, 1], and vector → a components vary from 2 to 0, linearly, throughout the iterations.

Hunting
The α, β, and δ type wolves have better awareness about the possible prey's position. Consequently, the initial three best solutions are saved. The other search agents should update their locations according to the position of the leading search agents [40] using Equations (55)-(61).

Attacking Prey (Exploitation)
When attacking the prey, the value of → a is reduced, which decreases the variation of → A. If |A| < 1, then, the next location of the search agent will be closer to the prey.

Search for Prey (Exploration)
The search is guided according to the α, β, and δ type grey wolves' positions. They go in different directions to search for prey, and gather again to attack it. This divergence is modeled using |A| > 1, which allows the GWO to search all over the space by forcing the search agent to get away from the prey. The → C vector is another constituent of the GWO that helps exploration. It contains random values between 0 and 2 inclusive. This parameter provides random weights for prey to emphasize (C ≥ 1) or deemphasize (C < 1) the effect of prey in determining the distance in Equation (51). Consequently, the GWO exhibits a random behavior during optimization to avoid local optima and promote exploration.
The GWO intentionally requires → C to provide random values to accentuate exploitation/exploration during initial and final iterations. This helps if there is a stagnation of the local optima. C is not linearly decreased in comparison to A.

Hybrid IC-GWA Optimization Approach
Hybrid algorithms are created to increase the performance of an optimization algorithm. They combine the advantages of two or more algorithms. The HIC-GWA is a combination of two evolutionary algorithms where the GWO is used to enhance the exploration ability of the ICA as shown in Figure 1. to accentuate exploitation/exploration during initial and final iterations. This helps if there is a stagnation of the local optima. C is not linearly decreased in comparison to A.

Hybrid IC-GWA Optimization Approach
Hybrid algorithms are created to increase the performance of an optimization algorithm. They combine the advantages of two or more algorithms. The HIC-GWA is a combination of two evolutionary algorithms where the GWO is used to enhance the exploration ability of the ICA as shown in Figure 1. In this proposed approach, ICA is initialized first to solve the OPF optimization problem. The assimilation and revolution of colonies, imperialist competition, elimination, and uniting empires are performed. The best solution of ICA is calculated as an initial condition of the GWA. The solution of the GWA is saved as the best value if it is less than the ICA's solution. The simulation continue until the stop condition is satisfied. The converged answer is achieved after termination of the algorithm.
The following steps show how to use the HIC-GWA to solve the OPF problem: i.
The power system data is specified. The HIC-GWA parameters are determined. ii.
Initialize the countries randomly, calculate their costs, and use assimilation. iii.
Exchange positions between imperialist and colony if it has a lower cost. v.
Calculate the total cost of all empires. vii.
Use solution obtained by ICA as initial condition for GWA. x.
The lower solution between ICA and GWA is saved as best solution. xi.
Go to step (ii) if the stop condition is not satisfied, otherwise, finish simulation. In this proposed approach, ICA is initialized first to solve the OPF optimization problem. The assimilation and revolution of colonies, imperialist competition, elimination, and uniting empires are performed. The best solution of ICA is calculated as an initial condition of the GWA. The solution of the GWA is saved as the best value if it is less than the ICA's solution. The simulation continue until the stop condition is satisfied. The converged answer is achieved after termination of the algorithm.
The following steps show how to use the HIC-GWA to solve the OPF problem: i.
The power system data is specified. The HIC-GWA parameters are determined. ii.
Initialize the countries randomly, calculate their costs, and use assimilation. iii.
Exchange positions between imperialist and colony if it has a lower cost. v.
Calculate the total cost of all empires. vii.
Use solution obtained by ICA as initial condition for GWA. x.
The lower solution between ICA and GWA is saved as best solution. xi.
Go to step (ii) if the stop condition is not satisfied, otherwise, finish simulation.

Simulation Results
The HIC-GWA has been applied on the IEEE 30 and 118 bus power systems to solve 12 different cases of OPF problems. The maximum number of iterations is 500 for IEEE 118 bus power system, and 100 for the IEEE 30 bus power systems. Power systems parameters are given in Table 1. The setting of the proposed HIC-GWA approach can be found in Table 2. MATLAB 8.3 (R2014a) has been used to implement simulations on a personal computer with i7 CPU 3.0 GHz 8.0 GB RAM [45,46]. The initial population is represented by N country . Each population contains one vector with N var components, including bus voltage and active power of the power plant, transformer tap changers, and shunt power injection compensator. The parameter N var , given in Table 1 is different for each case.
Solutions using the proposed approach will be compared with recently published OPF solutions using different optimization methods and objective functions shown in Table 3.

IEEE 30 Bus Test System
This power test system is used to exhibit the efficiency of the HIC-GWA. The details for busses and line data are shown in [43]. The system active and reactive power are 283.4 MW and 126.2 MVAR. The first five case studies have been used to solve simple objective OPF problems.

Case 1: Fuel Cost
This first single objective function considers minimizing the total fuel cost of power generation. It is modeled by the quadratic cost curve given in Equation (14). Simulation results, illustrated in Table 4, show that the fuel cost using the HIC-GWA is 798.20 ($/h). Compared with solutions from state-of-the-art existing optimization approaches in Table 5, the proposed HIC-GWA has significantly reduced the total fuel cost. The convergence curve of the total cost ($/h) for case 1 is shown in Figure 2. Note that it converged in less than 30 iterations. The first five case studies have been used to solve simple objective OPF problems.

Case 1: Fuel Cost
This first single objective function considers minimizing the total fuel cost of power generation. It is modeled by the quadratic cost curve given in Equation (14). Simulation results, illustrated in Table 4, show that the fuel cost using the HIC-GWA is 798.20 ($/h). Compared with solutions from state-of-the-art existing optimization approaches in Table 5, the proposed HIC-GWA has significantly reduced the total fuel cost. The convergence curve of the total cost ($/h) for case 1 is shown in Figure 2. Note that it converged in less than 30 iterations.

Case 2: Piecewise quadratic fuel cost
Thermal generators produce electricity by burning fuels such as coal, petroleum, or natural gas. The model for the fuel cost curve is given by Equation (15). Simulation results, illustrated in Table 4,

Case 2: Piecewise quadratic fuel cost
Thermal generators produce electricity by burning fuels such as coal, petroleum, or natural gas. The model for the fuel cost curve is given by Equation (15). Simulation results, illustrated in Table 4, show that the fuel cost using the proposed approach is 645.85 ($/h). Compared with existing optimization methods in Table 6, HIC-GWA has significantly reduced the total fuel cost. In cases 1 and 2, the proposed metaheuristic has a better convergence than recently published optimization methods.

Case 3: Piecewise quadratic fuel cost with valve point loading
The valve point loading effect is included in the cost function of Equation (17). Simulation results, illustrated in Table 4, show that the fuel cost using HIC-GWA is 902.25 ($/h). Compared with existing optimization methods in Table 7, HIC-GWA has significantly reduced the fuel cost in this case.

Case 4: Emission
The objective, in this case, is to reduce the emission level of important air pollutants like NOx and SOx, using the emission function described in Equation (18). Simulation results, illustrated in Table 4, show that the emission using HIC-GWA is 0.2009 (ton/h). Compared with existing optimization methods in Table 8, HIC-GWA has significantly reduced the emission level.

Case 5: Active power loss
To reduce transmission lines active power loss, we use the objective function given in Equation (19). Simulation results, illustrated in Table 4, show that the power loss using HIC-GWA is 2.61 (MW). Compared with existing optimization methods in Table 9, HIC-GWA has significantly reduced the power loss. In cases 3, 4, and 5, the proposed metaheuristic showed a better exploration than recently published optimization methods that appear to be stuck at a local minimum.

Multiobjective OPF
In the next five cases, we used the proposed metaheuristics to find better solutions for multiobjective OPF problems. Table 10 summarizes the best solutions of the simulation results using the HIC-GWA approach for cases 6-10.

Case 6: Fuel cost and active power losses
Cases 1 and 5 have been combined to reduce the fuel cost and the active power losses using the multiobjective function given in Equation (20). Simulation results show that the fuel cost and power loss using HIC-GWA are 856.99 ($/h) and 4.04 (MW). Compared with MSA, MDE, MPSO, FPA, and MFO approaches in Table 11, HIC-GWA has significantly reduced the fuel cost and power loss.

Case 7: Fuel cost and voltage deviation
Voltage profile management is essential to ensure system security. Voltage profile improvement reduces the deviation of load bus voltage. A multiobjective function is presented in Equation (21) to reduce the voltage deviations and fuel cost simultaneously. Simulation results show that the fuel cost and voltage deviations using the proposed approach are 802.45 ($/h) and 0.10 (p.u), respectively. Compared with MSA, MDE, MPSO, FPA, and MFO approaches in Table 12, HIC-GWA has significantly reduced the fuel cost and voltage deviations.

Case 8: Fuel cost with voltage stability improvement
The L index describes the system stability by measuring the distance of the actual state of the system to the stability limit. We are using the objective function given in Equation (22) to reduce both fuel cost and voltage stability. Simulation results, illustrated in Table 13, show that the fuel cost and L index using the proposed approach are 797.80 ($/h) and 0.11 (p.u), respectively. Compared with MSA, MDE, MPSO, FPA, and MFO approaches in Table 13, HIC-GWA has significantly reduced the fuel cost and L index.

Case 9: Fuel cost with voltage stability improvement during contingency condition
We consider the previous case with disruption of line (2-6) to simulate N -1 contingency. Best solutions for the fuel cost and the L index using HIC-GWA are 802.00 ($/h) and 0.11 (p.u), respectively. Compared with MSA, MDE, MPSO, FPA, and MFO approaches illustrated in Table 14, HIC-GWA has significantly reduced the fuel cost and L index during contingency condition. In cases 6-10, the proposed metaheuristic showed a better exploration than recently published optimization methods that appear to be stuck at a local minimum.
The total cost convergence curve for case 10 is displayed in Figure 3. The HIC-GWA approach converged in less than 50 iterations. Convergence curves of the fuel cost, voltage deviation, power loss, and emission are shown in Figure 4.

The IEEE 118 Bus Power System
The IEEE 118 bus test system [44], has been used for the next two cases to confirm the effectiveness of the HIC-GWA approach. The active and reactive power demand are 4242 MW and 1439 MVAR.

Case 11: Fuel cost
The function modeled by the quadratic cost curve given in Equation (14) is considered to minimize the total fuel cost of power generation. The simulation results, illustrated in Table 16, show Convergence curves of the fuel cost, voltage deviation, power loss, and emission are shown in Figure 4. Convergence curves of the fuel cost, voltage deviation, power loss, and emission are shown in Figure 4.

The IEEE 118 Bus Power System
The IEEE 118 bus test system [44], has been used for the next two cases to confirm the effectiveness of the HIC-GWA approach. The active and reactive power demand are 4242 MW and 1439 MVAR.

The IEEE 118 Bus Power System
The IEEE 118 bus test system [44], has been used for the next two cases to confirm the effectiveness of the HIC-GWA approach. The active and reactive power demand are 4242 MW and 1439 MVAR.

Case 11: Fuel cost
The function modeled by the quadratic cost curve given in Equation (14) is considered to minimize the total fuel cost of power generation. The simulation results, illustrated in Table 16, show that the HIC-GWA has significantly reduced the fuel cost compared with MSA, MDE, MPSO, FPA, and MFO approaches. In this case, the proposed metaheuristic has a better convergence than recently published optimization methods.

Case 12: Fuel cost with renewable energy sources (Wind/PV)
The objective in this case is to use the HIC-GWA to minimize the fuel cost (F 1 ), wind cost (F 2 ), and PV cost (F 3 ) for a system that includes renewable sources like WE and PV. The conventional power plants 12, 31, 66, 72, and 100 are replaced by five wind power units, and the conventional power plants 34,36,46, and 62 are replaced by four PV units. The simulation results are illustrated in Table 17. The total cost convergence curve for case 12 is presented in Figure 5. The proposed HIC-GWA approach converged in less than 100 iterations. In this case, the proposed metaheuristic has a better convergence than recently published optimization methods.

Case 12: Fuel cost with renewable energy sources (Wind/PV)
The objective in this case is to use the HIC-GWA to minimize the fuel cost ( ), wind cost ( ), and PV cost ( ) for a system that includes renewable sources like WE and PV. The conventional power plants 12, 31, 66, 72, and 100 are replaced by five wind power units, and the conventional power plants 34, 36, 46, and 62 are replaced by four PV units. The simulation results are illustrated in Table 17. The total cost convergence curve for case 12 is presented in Figure 5. The proposed HIC-GWA approach converged in less than 100 iterations.

HIC-GWA Robustness Analysis
Robustness analysis, which is a non-empirical form of confirmation, is an indispensable procedure in studying complex phenomena. A sensitivity analysis for case studies 1 and 11 has been performed to evaluate the robustness of the considered metaheuristic. Each parameter of the HIC-

HIC-GWA Robustness Analysis
Robustness analysis, which is a non-empirical form of confirmation, is an indispensable procedure in studying complex phenomena. A sensitivity analysis for case studies 1 and 11 has been performed to evaluate the robustness of the considered metaheuristic. Each parameter of the HIC-GWA has been perturbed by changing the values up and down. Likewise, optimization parameters values have been changed also to check the global effect of parameter's variations on the solution of the OPF problem. The equivalent Pareto solutions are illustrated in Table 18. The deviation ratio between normal and perturbed solutions is calculated using the following equation: Small deviations affirm the robustness of the HIC-GWA to variation of parameters in solving OPF problems.
To confirm the robustness of the HIC-GWA, we compare best and worst fuel cost averages to recently published OPF optimization methods in Table 19. The proposed HIC-GWA has consistently better solutions over 30 trial runs.  Table 20 shows the convergence speed of the HIC-GWA compared to recently published optimization methods. With 14.34 (s), HIC-GWA is second fastest to MFO by one hundredth of a second.

Conclusions
A novel hybrid optimization method combining imperialist competitive and grey wolf algorithm, HIC-GWA, has been proposed, developed, and applied successfully to solve twelve different test cases of single and multiobjective OPF problems in two IEEE test power systems with a mixture of wind energy and photovoltaic units. The results show that this metaheuristic is found to be very effective for large-scale applications, due to fast convergence and very few chances to get stuck at local minima. Analysis of the obtained solutions, along with a comparative study with recently published OPF optimization algorithms, proved the validity, effectiveness, and robustness of the HIC-GWA in precisely providing a set of stable optimal solutions, computed under realistic conditions, for a hybrid power system. This is very important in managing modern power systems, which are incorporating an ever-increased number of alternative energy sources. The proposed metaheuristic outperformed current well known and powerful algorithms in the literature, which confirms its superiority and potential to find valid and accurate solutions for multiobjective optimization. Indeed, the proposed paradigm may be used as a tool to answer many specific features of large-scale complex systems in general, thereby motivating further studies.