Uncertain Scenario Based MicroGrid Optimization via Hybrid Levy Particle Swarm Variable Neighborhood Search Optimization (HL_PS_VNSO)

Within the MicroGrid environment, the Energy Resource Management (ERM) problem becomes highly complex due to the uncertainty related to the Renewable Generation (RG) such as Photovoltaic power generation (PV), Electric Vehicle (EV) trip with Grid to Vehicle (G2V) or Vehicle to Grid (V2G), Energy Market price and load demand with Demand Response (DR) programs. Each of these issues should be tackled while optimizing revenues and reducing the running costs of Virtual Power Player (VPP) that collects each of these types of energy resources from the MicroGrid. This article presents a new hybrid optimization algorithm called “Hybrid Levy Particle Swarm Variable Neighborhood Search Optimization” (HL_PS_VNSO) to solve the ERM problem. Its key aspect is the hybridization of the Particle Swarm Optimization (PSO) and the Variable Neighborhood Search Optimization (VNS) algorithm with the enhanced step length using Levy Flight. The effectiveness of the proposed approach is measured by a 25-bus MicroGrid with 500 uncertain scenarios of the aforementioned uncertainty. The results of HL_PS_VNSO are compared with the most advanced optimization algorithms. The findings show that HL_PS_VNSO’s results are superior for the Average Ranking Index (A.R.I) and Ranking Index (R.I). For effective comparative analysis of algorithms, the traditional statistical method called One-way ANOVA Tukey Analysis is used. The results from this analysis also prove the superiority of HL_PS_VNSO over the most advanced optimization algorithms.


I. INTRODUCTION
In the MicroGrid context, the growing emergence of uncertain Distributed Energy Resources (DERs) such as solar PV, Demand Response (DR) programs, Electric Vehicles with G2V or V2G feature, Energy Storage Systems (ESSs) and the uncertain market price of energy are challenging the functioning of distribution networks. As a consequence, energy aggregators or Virtual Power Player (VPP) must tackle uncertain sources in a real-world microgrid problem. To do this, VPP needs a powerful technique to cope with the increasing variety of uncertainties. When important actions have to The associate editor coordinating the review of this manuscript and approving it for publication was Canbing Li . be undertaken every day in advance to optimize profit by lowering operating costs, day-to-day energy planning is a great challenge in the management of energy sources. However, owing to the vast number of energy sources and their inherent uncertainties, the ERM problem transforms into a Mix-integer nonlinear problem (MINLP) [1], which makes it quite complicated.
In the last decade, many Evolutionary Computation Techniques (ECT) have been widely used to solve the ERM problem with different uncertainty. In [2] Bacterial Foraging (BF) used for the MicroGrid optimization with the consideration of the cost of emissions as well as the operation and maintenance. Energy management in hybrid energy systems is discussed using the Interior Search Algorithm (ISA) [3].
Some of the popular ECT like Simulated Annealing [4], Particle Swarm Optimization (PSO) [5], Crow Search Algorithm (CSA) [6], Genetic Algorithm (GA) [7]- [10] and Tabu search (TS) [8] are used for ERM problem with different uncertainty. Artificial Bee Colony suggested [11] to address the day-ahead ERM in MicroGrid by taking into account uncertainties related to RG, EVs trip, market price and load demand. The Firefly Algorithm (FA) proposed in [12], for the Economical scheduling with optimized battery sizing. In [13] optimal scheduling is done by Imperialist Competitive Algorithm (ICA) in MicroGrid environment with uncertainty related to RG and load demand. Home Energy Management problem with uncertainty related to electricity price and solar PV is solved by Natural Aggregation Algorithm (NAA) [14]. Microgrid optimization with uncertainty related to RE source and cogeneration is achieved using the Efficient Heuristic Algorithm for Scheduling Energy Production (CHASE) [15]. Recently, hybrid or modified algorithms are becoming popular because of its competence to solve real-world complex problems. Signaled PSO (SiPSO) [16], Which is the updated variant of PSO with a signaling mechanism to tackle the ERM problem in the MicroGrid environment. In [17] the basic Gravitational Search Algorithm (GSA) is updated by self-mutation to prevent getting stuck in the local optima to solve the problem of energy resource scheduling. Adaptive Modified Firefly Optimization Algorithm (AMFA) [18], a refinement of the Firefly Algorithm utilizing a self-adapting approach to accelerate global space exploration. Enhanced variants of PSO called Gaussian Mutated PSO (PSO-MUT) [19] and Application Specific Modified Particle Swarm Optimization (ASMPSO) [20] are used to address MicroGrid scheduling in different scenarios. Modified Bacterial Foraging Optimization (MBFO) [21] technique suggested minimizing the operating cost as well as emission simultaneously. Hybridization of Simulated Annealing (SA) with initialization using Ant Colony Optimization (ACO) is proposed to tackle the ERM issue with high penetration of EVs [22]. Multi-period Artificial Bee Colony (MABC) [23] with uncertainty prediction using a Markov chain is applied for energy management. In [24], the performance of Differential Search (DS) and Quantum Particle Swarm Optimization (QPSO) is compared with its hybrid called HDSA and HQPSO for the energy management problem. The comparison reveals that HDSA gives better result in terms of profit and operation time as compared to the aforementioned algorithms. The modified version of PSO called Guaranteed convergence Particle Swarm Optimization with Gaussian Mutation (GPSO-GM) [25] and Adaptive Modified Particle Swarm Optimization algorithm (AMPSO) [26] are presented to solve the ERM problem. The GECAD (Polytechnic of Porto) group organized the competition on ''Optimal scheduling of distributed energy resources'' [27]. In this competition, Variable Neighborhood Search algorithm (VNS), Modified Chaotic Biogeography-based Optimization (CBBO) with Random Sinusoidal Migration and Cross-Entropy with Evolutionary Particle Swarm Optimization (CEEPSO) secured the first, second and third rank respectively. Grey Wolf Optimization with Fuzzy Logic (FL-GWO) is used for ERM and optimal battery sizing in microgrid environment [28]. The optimal allocation of renewable generation and energy storage system for economic benefits is done by Grey Wolf Optimization (GWO) [29]. In [30], Lezama proposed the hybrid-adaptive Differential Evolution (HyDE) algorithm for the ERM with uncertainties and compared it with the various version of DE. In MicroGrid with load and RE uncertainty, artificial fish swarm algorithm (AFSA) [31] is used for ERM. The Modified PSO suggested for real-time ERM with optimal usage of battery [32].
Hybridization of GA and SA is suggested in [33] for the optimal scheduling of the sources of energy in the islanded and grid-connected modes of a MicroGrid. In [34] author developed the hybridization of VNS and DEEPSO with cross-entropy characteristic for the optimal scheduling in a MicroGrid. The GECAD group in collaboration with Delft University organized the competition [35] on ERM problem in a microgrid with 100 uncertain scenarios at IEEE-WCCI/CEC conference, 2018. In this competition, Variable Neighborhood Search Optimization-Differential Evolutionary Particle Swarm Optimization (VNS-DEEPSO), Enhanced Velocity Differential Evolutionary Particle Swarm Optimization (EVDEPSO) and PSO with Global Best Perturbation (PSO-GBP)/Chaotic Evolutionary Particle Swarm Optimization Algorithm (CEPSO) got the first, second and third rank respectively. Table 1 provides the summary of the literature survey of the ECT techniques used for ERM problem with consideration of different uncertainty related to renewable generation (RG), forecasting load demand, Electric Vehicles (EVs) trips and Electricity Market (EM) price. In TABLE 1 the marks '' and 'x' indicate the consideration and no consideration of uncertainty, respectively in a respective reference.
Thus, the above-mentioned literature survey reveals the fact that, despite the modification and hybridization of algorithms, the objective of near-optimal solutions remains a major challenge in the field of optimization. No free lunch theorem [36] of optimization also proves that it is very difficult to find an algorithm that consistently provides the near-optimal and robust solution for all types of complex non-linear problems. Recently, hybridization of evolutionary algorithms has become popular in resolving the ERM problem due to their ability to handle non-linearity and uncertainty.
In summary, a robust hybrid optimization technique is key to achieving the near-optimal solution in a complex ERM problem with the highly uncertain environment. The article, therefore, has the following contributions to obtain a near-optimal solution: 1. Developed an effective hybridization of the PSO and VNS algorithm with the enhanced step length using Levy Flight entitled as ''Hybrid Levy Particle Swarm Variable Neighborhood Search Optimization''

II. PROBLEM FORMULATION OF ENERGY RESOURCE MANAGEMENT (ERM) WITH UNCERTAINTY A. VIRTUAL POWER PLAYER (VPP) ENERGY MANAGEMENT
The VPP aims to increase the profit by minimizing the operating cost of energy resources and maximize the income by selling the spare energy into the electricity markets. The objective function is converted into minimization of function f and represented by equation (1). Figure 1 illustrates the Energy Resource Management in MicroGrid by VPP with an uncertain environment. Where, VPP can buy the energy from the dispatchable DGs, renewable energy sources (i.e. PV units) and prosumer through the DR program. VPP can also buy or sell the energy from or to the energy storage units, Electricity Markets (EMs) and Electric Vehicles (EVs). After accumulating the energy from all the resources, optimal scheduling of energy is done by the VPP through the optimization and sells the energy to different types of load (i.e. residential, commercial and industrial load), EVs and EMs. The cost of buying the energy from the different energy sources is considered as an operating cost of VPP and the cost  (1) where 'Day' is the day before the optimal scheduling is done by VPP. O.C Day+1 total represents the total operating cost (O.C) of the resources managed by the VPP in a day advance. Total income generated by VPP is represented by Income Day+1 total .
The minimum value of f (hopefully negative) is the profit of the VPP. If f is negative it is expected to have a profit otherwise, O.C is higher than Income. Thus, the profit is P = −f , where P is the profit. Nevertheless, for the goal in optimization terms is to obtain the minimum value of f in the metaheuristics form (2), as shown at the bottom of the next page. In equation (2), O.C consists of the cost of dispatchable distributed generation (DG), external power suppliers, discharge of energy storage systems and electric vehicles, penalization of non-supplied demand, penalization of DG units' generation curtailment and Demand Response (DR) programs.
The indices are represented by: DG is an index of dispatchable distributed generation units; EPS is an index of external power suppliers; PV is an index of photovoltaic power generation units; ESS is an index of energy storage systems; EV is an index of electric vehicles; LD is an index of loads; SN is an index of scenarios; t is an index of a time slot.
The parameters are described by: N DG is the total number of DG units; N EPS is the total number of external power suppliers; N PV is the total number of photovoltaic power generation units; N ESS is the total number of energy storage units; N EV is the total number of electric vehicles; N SN is the total number of considered scenarios; N LD is the total number of loads. T is the total number of time slots. C DG(DG,t) is the cost of generation of DG unit in slot t (m.u./kWh); C EPS(EPS,t) is the cost of generation of EPS in slot t (m.u./kWh); C PV(PV,t) is the cost of PV units in slot t (m.u./kWh); C Dis(ESS,t) is the discharge cost of ESS in slot t (m.u./kWh); C Dis(EV,t) is the discharge cost of EV in slot t (m.u./kWh); C GC(DG,t) is the generation curtailment cost of DG unit in slot t (m.u./kWh); C NSD(LD,t) is the non-supplied demand cost of load LD in slot t (m.u./kWh); C CutDR(LD,t) is the DR load curtailment cost in slot t (m.u./kWh) The variables are described by: The Operating cost (m.u) is respresented as O.C; P DG(DG,t) is the active power generation by DG unit in slot t (kW); P EPS(EPS,t) is the active power generation by EPS in slot t (kW); Active power generation by PV unit in slot t for scenario SN (kW) is represented by P PV(PV,t,SN ) ; P Dis(ESS,t,SN ) is the active power discharge by ESS in slot t for scenario SN (kW); The active power discharge by EV in slot t for scenario SN (kW) is represented by P Dis(EV,t,SN ) ; P GC(DG,t,SN ) is the active P NSD(LD,t,SN ) power generation curtailment by DG unit in slot t for scenario SN (kW); P CutDR(LD,t,SN ) is the active power non-supplied demand of load LD in slot t for scenario SN (kW); P CutDR(LD,t,SN ) is the active power curtailment of load LD by DR program in slot t for scenario SN (kW); π (SN) is the probability of occurrence of the scenario SN. VOLUME 8, 2020 The VPP can receive its income by buying and selling of energy in the electricity markets [39].
In equation (3), Income consists of the cost of buying and selling of energy in the electricity markets. EM is an index of energy markets. N EM is the number of electricity markets. MP (EM,t,SN ) is the energy market price in slot t for scenario SN (m.u./kWh).
The variables are described by: Income is the VPP income (m.u.); P Buy(EM,t) is the active power buy from electricity market EM in slot t (kW); P Sell(EM,t) is the active power sell to electricity market EM in slot t (kW).
For the objective function (1) network considers the constraints of the power balance equations, EVs, ESSs, DR program and electricity market. For more details about the constraints please refer [39].

B. MODELING OF UNCERTAINTY
For the day ahead ERM, VPP completely relies on the forecast of uncertain sources like load demand, PV generation, market price and EVs trip. Assumption of an accurate forecast could create catastrophic failure of MicroGrid if actual data is different from the forecast. To cope up with this problem, the Monte Carlo Simulation (MCS) technique is applied for the scenario generations and reduction by considering the past data and trend of uncertain sources [40]. In MCS technique, the probability distribution function (PDF) is used to produce scenarios using the equation (4). Figure 2 shows the graphical representation of the scenario generation over each time slot t. where, x error,SN (t) is a normal distribution function with mean zero and standard deviation σ , x predict (t) is predicted value of variable x at t time and x SN (t) is the final value of variable x for scenario SN at t time. Also, in contexts of static metrics, the scenario curtailment method is used to exclude low-probability scenarios and aggregate the high-probability scenarios. Specific information on this technique is presented in [40]. The MCS technique generates the 5000 distinct uncertain scenarios of load production, PV generation and market price with forecast error using the normal distribution method of 10%, 15% and 20% respectively. The tool presented in [41] is used for the generation of uncertainty related to EVs trips. For this problem, the MCS approach reduced the scenarios to 500, depending on the highest probability of occurrence of scenario.

C. FITNESS FUNCTION
The fitness function presented in equation (5) is the sum of the objective function equation (1) and the aggregate penalty for violation of network constraints.
where, X represents the solution vector. c i is the magnitude of i th constraint, N CON is the total constraints and ρ is the penalty factor. Fit(X ) is the value of fitness function (m.u) for the solution vector X . By taking into account the uncertain sources in the solution vector, the fitness function is modified based on the scenario created by the MCS defined by equation (6).
where, δ SN is the alteration of variable in scenario SN, Fit SN (X ) represents the fitness function for the scenario SN. The fitness function mean and standard deviation [37] for the scenarios SN can be calculated using the equation (7) and (8) respectively: STD_Fit SN (X ) where, MEAN _Fit SN (X ) and STD_Fit SN (X ) are the mean and the standard deviation of fitness function over the 500 scenarios. The total number of times, the function evaluation is calculated by equation (9), which depends on the considered population, iteration and scenario for the solution of a problem. (9) where, N p is the total number of particles; N Iteration is the total number of Iterations; FEs is the total number of function evaluations. Figure 3 illustrates the internal functioning of the fitness function. The fitness function receives the input data as an array such as test case information, specific parameters, initial solutions and the number of scenarios (500 scenarios are considered for the competition). The internal functioning of the fitness function randomly chooses the 500 high probability scenarios from the total available scenarios and evaluate the fitness values of it. As an outcome, we get the fitness value of fitness function for each selected scenario as well as the solution of it in a matrix form.
The layout of the solution is an integral part of metaheuristics to display a given solution. The solution format adopted for this problem follows the vector representation as shown in Figure 4. The solution is represented as a vector of six-group of variables named as a group (1) active power output of DGs, group (2) generator binaries, group (3) EVs charge/discharge, group (4) ESSs charge/discharge, group (5) demand response of load and group (6) electricity market price. In each 1-hour duration, the vector comprising a total of 142 separate variables, which are reproduced chronologically for 24 hours. Thus, the total size of the solution vector for the 24 hours is 24 * 142 = 3408 variables. In a solution vector, all variables are continuous variables, apart from group (2) generator binaries. Group (2), binary of generator, refers to binary variables used to denote that a generator is coupled ('1'' value) or detached ('0' value). Further information about the internal functioning of fitness and solution are available in [37].

III. HYBRID LEVY PARTICLE SWARM VARIABLE NEIGHBORHOOD SEARCH OPTIMIZATION (HL_PS_VNSO)
The HL_PS_VNSO is the hybrid version of VNS and PSO with levy flight step length. In which VNS is used for the effective initialization of the population to get the near optimal solution in minimum execution time.

A. VNS ALGORITHM
The VNS [42] approach is to extend the local hunt and allow systemic improvements in the neighborhood and to address complicated problems. In VNS, the first task comprises the setting of the VNS strategic parameter and then initializing the population using equation (10).
where, X p,D is the initial position of the p th particle for dimension 'D' of the solution vector, Where, X MAX D and X MIN D is the upper and lower bound of vector 'D'. After the initialization of the particles, identify the best particle based on the minimum value of the fitness function using equation (5). After finding the best particle, call the Levy PSO. In which the best particle of VNS considers as an initial position.

B. LEVY PSO
Levy PSO is the modified version of PSO with levy distribution without inertia term. The new velocity in Levy PSO defined by equation (11). Where, the first, second and third step named as Perception, Cooperation and Levy step respectively. The product of the sum of all the three steps and Acceleration Factor_1(A.F_1) gives the new velocity of particle 'p'. Where A.F1 is used to increase the exploration of search space and U (0.2, 0.2, 1) represents the uniform random number between (0.2,1) with an interval of 0.2. The detail information about the three steps are given below:

1) PERCEPTION
In this step, the personal best particle P best followed the current position of particle X p . The perception term helps in the exploitation of the local search area to obtain the sub-optimal solution.

2) COOPERATION
In this concept, the global best particle G best followed by the current particle position. The G best is the particle, which has the minimum value of fitness function among all tested particles. Here, W C is the cooperation weight obtain by the random number between (0,1). In equation (11), W * C is the mutation of cooperation weight obtain by equation (12). The cooperation term useful for the exploration of global search space to obtain the sub-optimal solution.

3) LEVY STEP
It is a random walk, the length of which is derived from the Levy distribution as described in equation (13). Where, 'u' and 'v' obtain from the normal distribution. The most species (e.g. swordfish and Silky sharks) and insects use Levy flights to hunt for food [43]. In HL_PS_VNSO algorithm, the function of levy step is to efficiently exploit and explore the search space to obtain the global solution. The behavior of Levy flights in 50 successive steps beginning at origin (0,0) is illustrated in Figure 5.  where, u = rand(0, 1) * Sigma, v = rand(0, 1) where, β called levy coefficient.
The new velocity V new p given by equation (11) is the optimal step length of particles to get the optimal solution. The size of the step length is neither too low nor large; else, the solution gets stuck into local minima. Thus, to design the method in such a manner that the optimum value and direction of the step size would be given to achieve the sub-optimal solution in a sufficient amount of time. If the present velocity obtained by equation (11) exceeds the boundary point, it should be enforced by equation (15).
Here, the Limit Factor (L.F) enhance and keep the step length within the boundary limit.
After finding the new velocity V new p , determine the new position of each particle 'p' using equation (16). Equation (16) is the updated version of the standard new position equation as defined in all the population-based algorithms by using Acceleration Factor_2 (A.F_2) to speeds up the particle displacement process and prevent it from falling into the local minima. Figure 6 shows the movement of HL_PS_VNSO algorithm, in which the current position X p mutated by the A. is the vector sum of the term perception, cooperation and levy step. Figure 7 illustrates the flowchart of HL_PS_VNSO algorithm, which apply for the ERM problem in a highly uncertain environment. The flowchart of the HL_PS_VNSO algorithm has the following steps: Step 1 (Initialization by VNS): Here, VNS used for the effective initialization of the considered particles, for that first task comprises of setting the VNS strategic parameter and then initializing the particle population using equation (10).
Step 2: Set iteration to zero and determine the fitness of each initialized particle using the fitness function defined in equation (5) and find the best particle based on the minimum values of fitness function and then update the iteration.
Step 3: Improve the position of each particle by VNS to obtain the best particle.
Step 4: If iteration becomes greater than or equal to two, then call Levy PSO otherwise go to step-3.
Step 5 (Levy PSO): -Set the Levy PSO parameters and then initialize position and velocity between upper and lower bound by randomization.
Step 6: Determine the fitness function of each initialized particle and the best particle provided by VNS, and then update the memory based on the minimum fitness value. Update the iteration count by one.
Step 7: Replicate each particle's actual position, then update the memory by the replicated and best existing particle.
Step 8: Randomly produce a number between [0, 1] and if its value exceeds the probability of local search, go to Step-9 for exploring the global search space, otherwise go to Step-13 for the exploitation of local search space.
Step 10: Using equation (12), adjust the weight of the replicate population to generate the new particles.
Step 11: Impose the limit, if the current and replicate particles violating the limits.
Step 12: Determine the fitness function of each existing and replicated particle to obtain the global best particle and then go to step-16.
Step 13 (Local Search):-Using local search, determine the new location of the existing particles.
Step 14: Impose the limit, if the current particles violating the limits.
Step 15: Determine the fitness function of each existing particle to obtain the global best particle.
Step 16: Upgrade the memory of HL_PS_VNSO by the global best particle.
Step 17: Update the iteration, and then check the cutoff limit of iteration.
Step 18: If the cutoff limit is crossed, exit from the HL_PS_VNSO algorithm or else proceed to step-7.

IV. TEST CASE AND RESULT ANALYSIS
In this test case, the network shown is a real distribution system of a housing area in Portugal. The system is comprised of 24 underground lines linked with the main grid through the transformer at bus 1. All evaluated algorithms are tested for the 25-bus MicroGrid network contains 500 scenarios of uncertain sources such as PV units, 36 EVs, 90 residential loads with demand response scheme, 2 ESSs and 2 EMs. Table 2 shows the test case data of the 25-bus network in VOLUME 8, 2020   [39] terms of variable type, cost, power range and the number of units. Detail information about the case study available in [37]. Figure 8 represents the diagram of a 25-bus Micro-Grid network.
The VPP's objective is to optimally utilize the accessible distributed energy resources of the 25-bus MicroGrid. For that, HL_PS_VNSO algorithm was applied to solve the ERM problem under uncertain condition. To verify the efficacy and robustness of algorithms, for all competing algorithms we find the 20 final solutions (one for each run or trial) for the ERM problem. In each run, a maximum 50,000 function evaluations are allowed.
To prove the effectiveness of the HL_PS_VNSO algorithm, the results of HL_PS_VNSO are compared with the competition participated algorithm such as Gauss Mapped Variable Neighborhood Particle Swarm Optimization (GM_VNPSO),

CUMDANCauchy-C1: Cellular Estimation Distribution
Algorithm, PSO-GBP, Cross-Entropy Variable Neighborhood Differential Evolutionary Particle Swarm Optimization (CE_VNDEPSO) and EVDEPSO. The results of above all listed algorithms have been extracted from the competition database [37].
The HL_PS_VNSO algorithm is tested on a 64-bit operating system Intel Core(TM)-i3 processor with 8GB of RAM operating on Windows 10. MATLAB R2016b is used to address the HL_PS_VNSO algorithm.

A. FINE-TUNING OF HL_PS_VNSO PARAMETERS
One of the most critical aspects of the design of any metaheuristic technique is the calibration of the strategic parameters. In the HL_PS_VNSO method, a variety of experiments have been carried out to evaluate the optimum values of the strategic parameters. In these experiments, fix the size of population and no. of iteration to evaluate the fitness function closer to limit of 50,000 evaluations as per equation (9), but the other strategic parameters, such as Acceleration Factor-1 (A.F1), Acceleration Factor-2 (A.F2), Limit Factor (L.F) and local search probability altered to their optimum range. The results of the experiments give the optimum value of the strategic parameters as shown in Table 3 based on the best value of the Average Ranking Index (A.R.I).
After tuning the parameters of HL_PS_VNSO algorithm, it is used to solve the day ahead ERM problem. The performance of HL_PS_VNSO is compared with all tested VOLUME 8, 2020 FIGURE 8. 25-bus network [45]. algorithms in terms of average fitness and standard deviation for each run as shown in Table 4. The average fitness is the mean value of the fitness function over the 50,000 function evaluations in each run as per the equation (7). The standard deviation is the deviation of the fitness function over the 50,000 function evaluations in each run as per the equation (8). Table 5 presents the Ranking Index (R.I) of all evaluated algorithms as per the equation (17). The R.I is the sum of average fitness and standard deviation over 50,000 function evaluations in each run as per given in Table 4. In the case of metaheuristic optimization methods, the initial populations are created randomly. As a result, all the methods obtain different optimal solutions in every run (trial). So, for comparing the performance of various metaheuristic methods, the simulations are carried out many times (e.g. 20 runs) and the Ranking Index (R.I) are found out for each run. Table 5 reveals that HL_PS_VNSO has the best values of R.I in almost all 20 runs except 10, 11, 15, 16 and 17 runs. In runs 10, 11, 15, 16 and 17 the proposed HL_PS_VNSO algorithm gives inferior solutions as compared to the other two methods (GM_VNPSO and CE_VNDEPSO). So in 75% runs HL_PS_VNSO algorithm gives the best value of R.I compared to all tested algorithms. Out of 20 runs, the HL_PS_VNSO algorithm gives the best R.I in run 14 of value 78.53. It means that for this run HL_PS_VNSO gives the best solution for the ERM problem out of 20 runs. The best solution of run 14 given by the HL_PS_VNSO algorithm is provided in Table 10. This solution consists of the optimal scheduling of the different energy resources within the constraints limit. Table 6 shows the Average Ranking Index (A.R.I) of all tested algorithms for 20 runs obtained using equation (18). Here, N Runs is the number of runs used for the calculation of A.R.I. Table 6 also provides the average fitness, standard deviation, best and worst fitness function over 20 runs for all tested algorithms.    As seen in Table 6, the HL_PS_VNSO algorithm obtained the lowest A.R.I of value 84.10 m.u., which is the sum of the average fitness and the standard deviation. The minimum value of A.R.I. means that, the operating cost is less or profit is high. The other tested algorithms like GM_VNPSO, CE_VNDEPSO, CUMDAN Cauchy-C1, EVDEPSO and PSO-GBP got the second, third, fourth, fifth and sixth rank respectively in terms of A.R.I.
It is obvious from the comparison that HL_PS_VNSO delivers the best outcomes compared to all tested techniques in terms of A.R.I. The best and worst operating cost of HL_PS_VNSO is 3.76 m.u. and 284.93 m.u. respectively. The worst value of the operation cost obtained by HL_PS_VNSO is also the lowest compare to all tested algorithms. The standard deviation given in Table 6 indicates the variability VOLUME 8, 2020 of the fitness function over 20 runs. Out of all tested algorithms, HL_PS_VNSO provides the lowest value of the standard deviation of 47.52 m.u. Which shows the robustness of HL_PS_VNSO algorithm. Table 7 provides the assessment of the algorithms in terms of no of iterations and mean execution time. The execution time is also an essential aspect of checking the robustness of any algorithm. Table 7 reveals that, in terms of mean execution time, HL_PS_VNSO is ranked fourth among all tested algorithms. Whereas other algorithms such as CUM-DAN Cauchy-C1, CE_VNDEPSO, PSO-GBP, GM_VNPSO and EVDEPSO are ranked first, second, third, fifth and sixth respectively, but the A.R.I of them are quite high relative to HL_PS_VNSO. Therefore, an algorithm which provides the optimal solution in the lowest mean execution time would gain the benefit of solving this kind of large-scale problem in a limited time. The HL_PS_VNSO has discovered the best solution in a reasonable time, which indicates the efficacy of the suggested strategy.
However, a comparison focused solely on R.I indicates a weak approach to evaluating the results. Apart from the reality that the R.I value of HL_PS_VNSO is lower than all tested methods, it is important to verify whether all the methods have a statistically significant difference from each other in terms of R.I. For that, use the statistical test named as Oneway ANOVA Tukey's Honestly Significant Difference (HSD) test.

B. STATISTICAL TEST
One-way ANOVA [46] is a statistical approach used to verify the R.I of all the algorithms evaluated for each run shows any noticeable variations. In this respect, a test of the hypothesis is used to confirm the performance. In this test, the degree of significance is set at 5% to verify the statistical differences between the all tested algorithms. The one-way ANOVA test gives a meaningful finding, if the value of 'P' given by test is below 0.05, then it may be assumed that there is ample statistical proof that any one or more algorithms are substantially different from the other in terms of R.I. The limitation of One-way ANOVA is that it can only show that the R.I of all algorithms has substantial variances, but it cannot provide the details on which one algorithm differs from the others. For this purpose, a pairwise comparison test, known as Tukey's Honestly Significant Difference (HSD) test [47], is conducted  to categorize the which algorithms are substantially different from HL_PS_VNSO.
The one-way ANOVA analysis provides a value of 'P' below 0.05 as seen in Table 8, suggesting a major difference between one or more algorithms. Then after using the Tukey HSD test for pairwise assessment. In this test first step is to find the critical value (Q critical ) from the studentized range distribution table [48] based on the a = 6 treatments (algorithms) and degrees of freedom DF = 114 as given in Table 8. The critical value obtained from the studentized range distribution table is 3.172. Then, calculate the Qstatistic as per equation (19) for all the pairwise comparison with HL_PS_VNSO algorithm. The values of Q i,j are given in Table 9.
where, i, j = 1, ..a, i = j. y i − y j is the difference between the A.R.I of the compared pair of algorithms. M = 62.6 is the Mean Square for within algorithms as given in Table 8. Table 9 reveals that for all pairwise treatments, the magnitude of Tukey HSD Q statistic Q i,j is higher than Q critical except for the pairwise treatment between HL_PS_VNSO and GM_VNPSO. This implies that the R.I of HL_PS_VNSO and GM_VNPSO are approximately the same but substantially different from the remaining algorithms in each run. Whereas in terms of A.R.I, the HL_PS_VNSO algorithm provides the best outcome relative to all tested algorithms. Also as seen in Table 7, HL_PS_VNSO offers the solution in lower mean execution time as compared to GM_VNPSO. So overall the HL_PS_VNSO provides the best solution for ERM problem in a highly uncertain environment as compared to all tested methods.

V. CONCLUSION
In the Microgrid context, the extensive usage of uncertain sources of energy such as PV, ESSs, DR schemes, V2G and Electricity Markets with various scenarios making the ERM challenge highly dynamic and high-dimensional. To solve this problem with maximum profit, the VPP must use the appropriate method.
This article suggested a hybrid version of VNS and PSO with levy flight step length and entitled as Hybrid Levy Particle Swarm Variable Neighborhood Search Optimization'' (HL_PS_VNSO). In the proposed method, the VNS is used for the effective initialization of the population to reach towards the near-optimal solution in minimum iterations and execution time. The use of the levy step in the PSO algorithm enhances the global exploration capability of the proposed algorithm.
The HL_PS_VNSO algorithm is applied for the ERM problem under highly uncertain environment. The robustness and efficacy of HL_PS_VNSO are compared with the most advanced optimization algorithms namely GM_VNPSO, CE_VNDEPSO, CUMDAN Cauchy-C1, EVDEPSO and PSO-GBP. Comparative evaluation reveals that HL_PS_VNSO provides the low value of R.I and A.R.I relative to the algorithms described above. This means that VPP acquires maximum profit by using the HL_PS_VNSO algorithm. Statistical evaluation by the one-way ANOVA Tukey HSD test further confirms its superiority by demonstrating that the R.I of HL_PS_VNSO is statistically substantially different from the above-listed algorithms except for the GM_VNPSO algorithm. But in terms of A.R.I and mean time of execution, HL_PS_VNSO offers better performance than GM_VNPSO. Eventually, the suggested HL_PS_VNSO approach produces optimal performance and demonstrates that it is capable of coping with practical problems containing a wide variety of uncertainty.
In the future, HL_PS_VNSO will be used to tackle an extremely dynamic MicroGrid network with large integration of unpredictable energy sources and a broad range of scenarios.

APPENDIX
The optimal solution of the HL_PS_VNSO algorithm is provided in Table 10.
His team's proposed optimization methods entitled EVDEPSO and IC-DEEPSO had secured second and fourth ranks in the IEEE PES and CIS sponsored optimization competitions on Evolutionary Computation in Uncertain Environments: A Smart Grid Application at IEEE-WCCI conference 2018, Rio De Janeiro, Brazil. Also, his team's proposed optimization method entitled HL_PS_VNSO secured the second rank in the same competition title name but with more complexity in the problem at the IEEE CEC Conference 2019, Wellington, New Zealand, and the GECCO Conference 2019, Prague, Czech Republic. His research interests include artificial intelligence, energy resource management demand response in smart grid, microgrid, and machine learning application in the power systems. His research interests include power system optimization, computational intelligence methods, restructured power systems, smart grid, and renewable integrations. He has won many IEEE and CIS sponsored Smart Grid Optimization Competitions. His team's proposed optimization methods entitled Levy DEEPSO and EE-CMAES had secured third and second ranks at the IEEE PES Optimization Competitions at the IEEE PES General Meetings, USA, in 2017 and 2018, respectively. His team's proposed optimization methods entitled EVDEPSO and IC-DEEPSO had secured second and fourth ranks in the IEEE PES and CIS sponsored optimization competitions on Evolutionary Computation in Uncertain Environments: A Smart Grid Application at the IEEE-WCCI Conference 2018, Rio De Janeiro, Brazil. Also, his team's proposed optimization method entitled HL_PS_VNSO and GM_VNPSO secured second and third ranks, respectively, at the IEEE CEC Conference, New Zealand, and the GECCO 2019, Czech Republic, in 2019. VOLUME 8, 2020