Cross Entropy Covariance Matrix Adaptation Evolution Strategy Solving the Bi-Level Bidding Optimization Problem in Local Energy Markets with Renewables

Abstract: The increased penetration of renewables in distribution power systems has motivated researchers to take significant interest in local energy transactions. The major goal of Local Energy Markets (LEM) is to promote the participation of small consumers in energy transactions and providing an opportunity for transactive energy systems. Such energy transactions in LEM are considered as a bi-level optimization problem in which all agents at upper and lower levels try to maximize their profits. But typical bi-level problem is very complex as it is inherently nonlinear, discontinued and strongly NP-hard. So, this article proposes the application of hybridized Cross Entropy Covariance Matrix Adaptation Evolution Strategy (CE-CMAES) to tackle such a complex bi-level problem of LEM. The proposed CE-CMAES secured the 1st rank in Testbed-2 entitled, “Bi-level optimization of end-users’ bidding strategies in local energy markets (LM)” at international competitions on Smart Grid Problems, held at GECCO 2020 and WCCI 2020. CE method is used for global exploration of search space and CMAES is used for local exploitation as its adaptive step-size mechanism prevents its premature convergence. A practical distribution system with renewable energy penetration is considered for simulation. The comparative analysis shows that the overall cost, mean fitness and Ranking Index (R.I) obtained from CE-CMAES are superior to those obtained from the state-of-the-art participated algorithms. Wilcoxon Signed Rank Statistical test also proves that CE-CMAES is statistically different from the tested algorithms.


Introduction
The smart grid era has created lots of new challenges and opportunities that must be explored and implemented to optimize its performance [1]. The large penetration of distributed renewable energy sources has created significant opportunity in energy transactions at the local lower level [2]. In this paradigm, local energy markets (LEM) facilitates small agents to trade renewable energy at the local level and contribute to the reduction of environmental pollutants. Thus, LEM empowers them to take a large part in the energy community to implement fully transactive energy systems [3]. The energy transactions in LEM results in to a bi-level optimization problem in which all agents try to maximize their profits by optimizing their bidding strategies. In this regard, [4] proposed bi-level problem to optimize offering prices of energy producers and subsequently it was converted into a mixed integer linear programming. But demand side bidding was ignored and not optimized.Proposed optimal bidding strategy of transactive agents in which two subproblems were considered to reduce no. of variables and nonlinear constraints were re-placed by linear inequalities to reduce the computational burden [5]. But the effect of renewables was ignored in the problem formulation, which is important for modern distribution system. Authors of [6] proposed a two-stage robust model to optimize bidding in day-ahead and real time markets in virtual power plant environment. The optimization problem was formulated using linear programming. But since the optimization problem was linearized, the quality of the optimal solutions was also compromised. Authors of [7] proposed hybrid optimal bidding strategy for a microgrid in day-ahead market under the influence of intermittent distributed generation and price responsive loads. Stochastic models were used to model day-ahead market price and uncertain output of intermittent DG while MILP was used to limit power imbalance in real time market considering uncertainty in real time market price. The hourly bidding curves obtained by such approach were less accurate as the original nonlinear problem was considered as a linear. [8] Proposed a robust based bidding and offering strategy for a price-maker energy storage facility participating in a day-ahead market. It proposed max-min MILP approach to present a participation strategy to manage the risk of uncertainty associated with forecasted hourly generation and demand price quota curves. Some authors have explored metaheuristic optimization methods to solve bi-level problem in LEM. An example of such application can be found in [9], where Particle Swarm Optimization (PSO) was used to solve bi-level problem in which power producers' biddings were optimized at upper level in monthly contract and balancing markets. In this, a basic PSO is implemented whose local exploitation search ability is poor and thus suffers with pre-mature convergence. Authors of [10] used basic Genetic Algorithm (GA) to optimize the bidding of price maker retailers with probability-based estimation of market price. But the GA is quite computationally expensive. In [11] Ant colony optimization was proposed to optimize biddings of all market participants in LEM.
The literature survey reveals that the bidding optimization problems of the market participants are generally formulated as a bi-level optimization problem. In the majority of the surveyed papers, mixed integer linear programming method was used. In this, the objective functions/constraints are linearized and as a result the accuracy of the optimal solutions is decreased. But typically, any bi-level optimization problem is inherently highly nonlinear [12], non-convex, non-differentiable, discontinued and strongly NP-hard [13], because the lower-level problem has nonlinear inducible region and it may have multiple optimal solutions as shown in Figure. 1. Also, upper and lower level problems are mutually dependent. Remaining surveyed papers have used basic meta-heuristic methods to solve bi-level problems. Such methods are either vulnerable to pre-mature convergence or computationally expensive.
So, to get rid of aforementioned problems related to bi-level optimization, this article provides the following contributions to the available literature.
• A multi-period bi-level optimization problem is formulated in which competitive agents in the upper-level problem try to maximize their profits, modifying and depending on the price determined in lower-level problem (i.e. the clearing price in the LEM), hence resulting in a strong interdependence of their decisions. • The applications of 1st rank winner optimization algorithm entitled "CE-CMAES" in Test Bed 2: "bi-level optimization of end-users' bidding strategies in local energy markets" at international competitions entitled, "Evolutionary Computation in Uncertain Environments: A Smart Grid Application" organized at Genetic and Evolutionary Computation Conference (GECCO 2020) and IEEE world Congress on Computational Intelligence (WCCI 2020) have been proposed to solve bi-level optimal bidding in local energy markets with renewable energy sources. • Practical distribution power system with renewable energy sources has been used to test the effectiveness of the proposed algorithm. • The comparative analysis shows that the profits of all agents obtained from CE-CMAES are better than those obtained from the state-of-the-art algorithms in a dayahead local energy market.
The remaining paper is organized in six sections. Following this introduction, the formulation of bi-level optimization problem is discussed in section 2. Section 3 discusses the CE-CMAES algorithm in detail. In section 4, the comparative analysis of the proposed method with the state-of-the-art optimization methods is carried out on a practical distribution system. Eventually, section 5 outlines main conclusions and prospective future research work. In section 6 references are given.

Optimization Problem Formulation
A day-ahead Local Energy Market (LEM) bidding optimization problem is proposed in [14] in which agents submit offers and bids to maximize their profits and minimize their costs respectively. The agents considered in LEM are of the 3 types, namely prosumers (consumers with power generation capabilities), small producers and consumers. Also, all agents have access to the main grid, which works as a back-up power system. So, all agents can trade energy in LEM with prices between the feed-in tariff ( ) F c and Grid electricity tariff ( ) G c . It is assumed that FG cc  . So, it is assumed that direct energy transactions with the main grid are less beneficial to the agents than in the LEM. Figure. 2 shows the LEM scenario.

Bi-level optimization problem formulation
The LEM bidding optimization problem is basically a bi-level optimization problem in which the upper-level problem and lower-level problem correspond to the maximization of agents' profits and maximization of energy transacted respectively in LEM [14]. So, the solution of the lower-level affects the upper-level by modifying the profits of all agents. where each agent j wants to maximize its profit and each agent i wants to minimize its costs.

Upper Level Problem
The optimization problem for all producer agents can be formulated as Equation (1) in which they try to maximize their profits regarding their production marginal cost. ,, Where, j P is the profit of agent j , cp is the LEM clearing price (equal for buyers and sellers), , ji x is the energy sold by agent Esell is the energy sold by agent j to grid, and mc j cG is the marginal cost associated to j . mc c is taken as zero for PV generation. Let, is the marginal cost associated to a Combined Heat and Power (CHP) generator and it is defined as a monotonically decreasing function [15] as given in Equation (2).
Where, CHP b is a constant cost factor of the CHP generation unit and j G is the energy produced by the CHP unit.
The optimization problem for all consumer agents can be formulated as Equation (3) in which they try to minimize their costs. ,,

Lower Level Problem
It is evident that the agents' profits and costs are dependent upon the LEM clearing price, which is determined in the lower-level problem and depends upon the bidding process. The lower-level problem is modelled as an asymmetric pool market in which offers and bids are allocated using a merit order mechanism to determine the supply and demand curves for energy [6]. It determines the clearing price ( ) cp and energy transacted between the agents. A simple deterministic mechanism is used to determine the amount of energy traded to the grid (i.e., variables Where, marginal, j P are the profits that agent j could make considering the marginal cost of producing its maximum capacity and selling the non-traded energy in the LEM to the grid. j P are the regular profits that the agent would make by selling only into the LEM. As a result, it is guaranteed that producers only sell energy to the grid if they can actually get some profits with their remaining capacity. In case of PV generation where the marginal cost is almost zero, the agents will always sell its remaining energy to the grid as per this rule. Consumer i decides the energy bought from the grid as per Equation (5 x is the energy bought by consumer i from producer j in LEM. So, using Equation (4) and Equation (5), the second terms of Equation (1) and Equation (3) respectively could be found out.

Formulation of objective function
In Equation (1) and Equation (3), the clearing price is the equilibrium price at which supply and demand become same. The bid of each agent influences the market price, and thus a strong inter-dependency exists between profits and costs of agents. The aim of optimization is to maximize the overall average profits of the power system and to provide optimal solutions that distribute the earnings among all agents. Thus, the objective function is formulated as Equation (6) (6) Where profit ( ) j P and cost ( ) i C are conflicting objectives since all agents try to achieve the best results for their own.

Representation of solutions and fitness function
As per Equation (6), Bi-level optimization tries to maximize profit of all agents by executing optimal bidding of agents in the LEM. Let, is a set which consists of all agents (producers and consumers). The aim is to find the best tuple ( ) , kk q p k K  representing the quantity of active power to bid and optimal price respectively in the LEM for each agent. The bidding is done for all tT  periods of optimization. For day-ahead market 24 T = periods are taken. So, the vector of optimization variables x is represented as , kt q is the quantity and , kt p is the price the th k agent will submit to the LEM at time period t . A sign convention is used to define consumer and producer type. A positive quantity represents a bid (i.e., buying in the LEM), and a negative quantity represents an offer (i.e., selling in the LEM). Thus, the action of an agent can be controlled by defining the bounds of variables in which consumer agents can submit bids in the LEM within the bounds   max 0, L (i.e., between 0 and their maximum consumption limit), while producer agents can submit their offers within the bounds max, 0 PC  −  (i.e., between 0 and their maximum power production capacity). The minimum and maximum boundary values for prices are in the range , FG cc   and same for all agents. Figure. 3 shows the structure of population of particles (solutions). Such particles are evaluated in the fitness function that returns the mean average profit of all agents plus the standard deviation as per Equation (7).
Where ( ) mean profits and ( ) std profits are functions that compute the average and standard deviation of the profits that all agents obtained considering the bids/offers encoded in the particle. The negative sign in the first term is used to transform the profit maximization problem into a minimization problem. The less the value in Equation (7), the better the mean profits achieved by all agents. Equation (7) is optimized using the proposed CE-CMAES method, which is described, in the following section.

CE-CMAES Method
CE-CMAES is a sequential hybridization of Cross Entropy (CE) proposed by Rubinstein [16] and CMAES was proposed by Nikolaus Hansen [17]. It has been applied to solve dynamic optimal power flow [18]. Both are population-based iterative meta-heuristic optimization methods similar to Particle Swarm Optimization [19] and Differential Evolutionary methods [20]. However, unlike conventional meta-heuristic methods which directly act upon prospective solutions (particles) to obtain optimal solutions, these methods generate "distribution" of prospective solutions and update it by updating the parameters of the distribution to obtain the optimal solutions. CE method was used for global exploration of the search space in the initial part (50% iterations). As it has small amount of parameters to be tuned, it can be easily implemented. CMAES method was used for local exploitation of search space for the remaining 50% iterations. It updates the mean of population such that the likelihood of elite particles near the sub-optimal solutions is increased. Covariance matrix is updated such that the likelihood of good points with better fitness and search steps is increased. Also, the evolution path provides a unique adaptive step size which prevents its premature convergence and thus enhances its local exploitation capability. The following section explains the theory of CE method.

Cross Entropy (CE) Method
The bids of all agents are optimized to minimize Equation (7). Firstly, a set of particles are randomly generated which obey the probability distribution function (pdf) ( ) where  is the vector of parameters to be optimized. Generally, ( ) .; f  is the Gaussian distribution parameterized by its mean m and variance 2  . i.e. f  . The following section shows the step-by-step procedure of CE-CMAES.
Step: 1 Initialization of mean and standard deviation ,c   for bid prices.) Mean and standard deviation of population are initialized as per Equation (8) and Equation (9) respectively.
The conventional method to generate ( Step: 2 Generation of population of particles Generate the population of particles from the sampling distribution gg Nm  as per Equation (10).
Where, a NP is the total no. of particles (Population Size), ( 1) is the mean value of the search distribution for each dimension at generation g , ( randn is a normally distributed random variable with parameters Step: 3 Adjustment of search distribution limit violations The maximum and minimum boundary limits of all particles are checked and set as per Equation (11).
Step: 4 Execution of Lower-level problem Particles' bids and offers are used to determine market clearing price ( ) cp through merit order mechanism. Obtain energy transacted between the agents using Equation (4) and Equation (5).

Step: 5 Execution of Upper-Level Problem
Calculate the profit/cost of all agents using Equation (1) and Equation (3).
Step: 6 Fitness Function Evaluations and ranking The fitness of all particles are evaluated by using Equation (7) which calculate the mean average profit of all agents and its standard deviation. Sort (rank) all fitness values in ascending order as given in Equation (12).
Where, () i fx is the fitness of th i particle, 1 x is a Global Best ( ) Best G particle having the minimum fitness among all particles. Consider top best 20% particles as "elite" particles.

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 8 December 2021
Step: 7 Updating the mean and standard deviation of elite particles The mean ( ) of the selected elite particles is found by Equation (13). x + is the th i best particle among whole population at ( NP denotes the index of the th i rank particle. Standard deviation ( ) ( ) 1 g   + of elite particles is found by Equation (14). As elite particles are near to sub-optimal solutions, more smoothing (weightage) is applied to their mean value as compared to mean of whole population as per Equation (15). Similarly, the standard deviation of elite particles should be modified to a small extent as they lie in the vicinity of sub-optimal solutions. So, less smoothing is applied to them as compared to standard deviation of all particles which require more exploration and thus more smoothing as given in Equation (16).
Step: 9 Increment of generation count Set :1 gg =+ Go to step #2 and repeat the steps 2-9 until 50% of max g are completed. So, CE method is used for global exploration. Subsequently, CMAES method is applied for local exploitation, which is discussed in the following section.

CMAES method
Step: 10 Initialization of mean of each dimension Global Best ( ) Best G particle obtained from CE method is considered as the mean value of corresponding dimension of CMAES as shown in Equation (17). As a result, the new search distributions (individuals) will be generated near to the sub-optimal solutions. So, the local exploitation of the search space becomes more effective as compared to with traditionally initialized mean value of each dimension.
is the covariance matrix at generation g .
Step: 12 Adjustment of search distribution limit violations Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 8 December 2021 The maximum and minimum boundary limits of all particles are checked and set as per Equation (11).
Step: 13 Execution of Lower-level problem Particles' bids and offers are used to determine market clearing price ( ) cp through merit order mechanism. Obtain energy transacted between the agents using Equation (4) and Equation (5).
Step: 14 Execution of Upper-Level Problem Calculate the profit/cost of all agents using Equation (1) and Equation (3) Step: 15 Fitness Function Evaluations and ranking The fitness of all particles are evaluated by using Equation (7) which calculate the profits/costs of all agents and the standard deviation between the profits/costs of agents. Sort (rank) all fitness values in ascending order and find Best G and top 50% elite particles ( )  which have better fitness than remaining 50% particles.
CMAES is an iterative method in which ( ) ( ) 11 , gg mC ++ and ( ) 1 g  + are updated in every generation to obtain better solutions.
Step: 16 Updating the Mean ( ) The updated mean of the search distribution is a weighted average of selected  elite particles from the whole population as per Equation (19).
x + is the th i best particle among whole population at ( 1) g + generation. The index : NP denotes the index of the th i rank particle.
1 w is the weightage applied to Best G particle. As the highest weightage is applied to Best G particle, the new individuals will be pulled towards it. Thus, update of mean maximizes the likelihood of elite particles near the sub-optimal solutions.
Step:17 Updating the Covariance Matrix ( ) ( ) 1 g C + A Covariance Matrix represents pairwise dependencies between the decision variables and it determines the shape of the distribution ellipsoid. It is updated such that the likelihood of good points with better fitness and search steps is increased. A covariance matrix of next generation is produced from the current covariance matrix as per Equation (20).
where, C  is a learning rate for updating the covariance matrix, 01 C   if 0 C  = : No learning takes place, if 1 C  = : No prior information is retained To enhance this step-size, the "evolution path" is also used. It is the sum of consecutive steps. To construct the evolution path ( ) 1 g c p + , the exponential smoothing is applied to assign more weightage to recent generations as per Equation (21). By using Equation (21), Rank-one update of the covariance matrix ( ) g C is obtained as per Equation (22). Where, The combination of (20) and (22)

 + The
Step-size is controlled by an evolution path, which is a sum of successive steps independent and uncorrelated. The length of evolution path is decided by comparing its length with its expected value. If selection biases the evolution path to be longer than expected, then step-size is increased and vice versa. Ideally, selection does not bias the length of the evolution path and the length equals its expected length. A conjugate evolu-

Start
Step

1: Initialization of Cross Entropy(CE) Method: Set no. of particles, Dimension of the problem, elite particles, Smoothing parameters, Initialize mean and standard deviation of each dimension using (8) and (9) respectively, generation (g)=0
Steps 2 & 3: Generate the population of particles using (10) and adjust their limit violations using (11)

Is g<0.5*gmax?
End: Display Optimal prices, quantities, and costs/profits of all agents

YES
Step 4:Execute the Lower level problem: Allocate offers and bids using a merit order mechanism to determine the supply and demand curves. Calculate the market clearing price. Calculate energy transacted between the agents using (4) and (5) Step 5: Execute the Upper level problem: Calculate the profits and costs of all agents using (1) and (3) respectively Step 6: Fitness function evaluations and ranking: Calculate average and standard deviation of profits of all agents by evaluating the fitness of all particles using (7). Rank all fitness values in ascending order using (12) Step 7: Obtain the Global Best particle. Find mean and standard deviation of elite particles using (13) and (14) respectively Step 8: Apply smoothing to mean and standard deviation of elite particles and whole population using (15) and (16) respectively Step 9: Increment: g=g+1

Cross Entropy method for Global Exploration
Step 10: Apply CMAES :

Global Best particle of CE is considered as the mean value of each dimension using (17), initialize Covariance matrix, and Initialize various parameters of CMAES
Steps 11 & 12: Generate the population of particles using (18) and adjust their limit violations using (11) Step 13: Execute the Lower level problem: Allocate offers and bids using a merit order mechanism to determine the supply and demand curves. Calculate the market clearing price. Calculate energy transacted between the agents using (4) and (5)

Step 14: Execute the Upper level problem: Calculate the profit and costs of all agents using (1) and (3) respectively
Step 15: Fitness function evaluations and ranking: Calculate average and standard deviation of profits of all agents by evaluating the fitness of all particles using (7). Rank all fitness values in ascending order using (12) Step 16: Update the weighted mean of whole population using (19) and assign highest weight to Global Best particle Step 17: Update the evolution path of Covariance matrix and Covariance matrix using (21) and (23) respectively Step 18: Update the evolution path of Step-size and Step-size using (24) and (26) respectively Is g<=gmax?
Step   The step-size (standard deviation) is adaptively changed which ultimately prevents premature convergence of CMAES [21] and improves local exploitation of solution space.
Step: 19 Increment of generation count and Convergence check Set :1 gg =+ Go to step #11 and execute steps 11-19 repeatedly until maximum number of generations are exhausted. The output of CE-CMAES are optimized prices, quantities (powers), and costs/profits of all agents. The flowchart of the Bi-level bidding optimization using CE-CMAES is given in Figure 4. The MATLAB™ codes of participated algorithms are available at http://www.gecad.isep.ipp.pt/ERM-competitions/2020-2/ [22].

Case Study and Result Analysis
In this section, the case study and the application of a set of algorithms to solve the bi-level optimal bidding problem in LEM are presented. We use the case study from the competition that considers nine agents, 3 of which are consumers, 3 prosumers (i.e., consumers with PV generation capabilities), and 3 are small generators (i.e., CHP generators). The case study data generated from the reference power profiles of residential houses and PV systems were developed using the open datasets available on the PES ISS website [23]. Figure 5 shows the base profiles and range of power of three normal houses and a PV power profile used to generate the data. To do this, a randomized function with a uniform distribution around 20% of those profiles was applied. We also included generator agents represented as small CHPs generators with a maximum generation power of 2 kW and a marginal cost determined using Equation (2) with a factor of 0.18 chp b = EUR/kWh [15]. At the end, feed-in and grid tariffs are set to 0.12 F c = and 0.28 G c = as in [15].
The purpose of the bi-level optimization is that all agents maximize their profits. With this purpose, the CE-CMAES algorithm was applied to solve the problem. To prove its effectiveness, the results of CE-CMAES were compared with 11 algorithms participating in the 2020 competition on "Evolutionary Computation in the Energy Domain: Smart Grid Applications" [22]. The 11 algorithms include Recursive Differential Grouping 3-Differential Evolutionary Particle Swarm Optimization (RDG3-DEEPSO), Ensembled method of CBBO, Cauchy and DEEPSO algorithm, Enhanced Hybrid Levy Particle Swarm Variable Neighborhood Search Optimization (EHL_PS_VNSO), CUMDANCau-chy++: a Cellular EDA [24], HFEABC, Differential Evolutionary and Estimation of Distribution Algorithm (DEEDA), Differential Evolution and Teaching Learning Based Optimization (DE-TLBO), GASAPSO, AJSO, Hybrid Adaptive Differential Evolution with Decay Function (HyDE-DF) and Particle Swarm Optimization with Global Best Perturbation (PSO-GBP). To check the effectiveness and robustness of algorithms, we considered the 20 final solutions (one for each run) of the problem for each of the competing algorithms. A limit of 50,000 function evaluations were considered in each run. The results of all the above-mentioned algorithms were collected from the competition database [22].

Fine-Tuning of CE-CMAES Parameters
One of the most critical aspects of the design of any metaheuristic technique is the calibration of the strategic parameters. The CE-CMAES method has only two smoothing parameters α and β to be tuned. Thus, a variety of experiments had been carried out to evaluate their optimum values. In these experiments, we fixed the population size 100 a NP = , No. of Iteration (gmax)=500 and changed the smoothing parameters α and β within their feasible range of (0.7,1) and (0,0.3) respectively. The results of the experiments gave the optimum values of the smoothing parameters α and β, which were 0.9 and 0.1 respectively based on the best value of the mean fitness as shown in Table 1.

Performance Comparison and Results
After tuning the CE-CMAES smoothing parameters, the algorithm was applied to solve the bi-level optimization problem. To measure the CE-CMAES algorithm's performance, it was compared with contemporary state-of-the-art optimization algorithms in terms of fitness value for each run, as shown in Table 2. The fitness value displayed in Table 2 is the average value of the fitness function over the 50,000 function evaluations of each run. runs. It reveals the robustness of the CE-CMAES algorithm relative to other tested algorithms. Table 2 also reveals that the 11th run of CE-CMAES offers the best value of fitness of 2.0988 EUR out of 20 runs. It means that, in this run, CE-CMAES provides the best solution to the bi-level optimization problem. The optimal bids/offers and prices of all agents obtained by CE-CMAES in the eleventh run is given in Appendix A (Table A1). All algorithms were run 20 times, and the results were recorded to obtain the mean value of those 20 runs. Table 3 gives the total overall costs of the system, the costs/profits of the group of agents (i.e., producers, prosumers and consumers), the mean fitness, the standard deviation, R.I. and the time taken by each algorithm to execute the 20 runs. The mean fitness is the average fitness over 20 runs (trials). The Ranking Index (R.I) is the sum of the mean fitness and the standard deviation over 20 runs. Table 3 also includes the overall costs/profits that agents can attain in the absence of a LEM (used as a baseline). The CE-CMAES algorithm achieved the first rank in terms of R.I with the lowest value of 2.1148. Also, EHL_PS_VNSO and CUMDANCauchy++: a Cellular EDA achieved the second and third rank with R.I of 2.1491 and 2.1506, respectively. Notice also that CE-CMAES achieves the best overall costs of 3.1840 EUR for the system compared to all tested algorithms. In terms of execution time, CE-CMAES takes an average of 17.0210 minutes to solve the problem, which is lower than the EHL_PS_VNSO (25.4820 minutes). It proves the advantage of hybridization of CMAES with CE, as the best solution obtained by CE is used for the initialization of CMAES. As a result, CMAES is near to the sub-optimal solutions. Therefore, CE-CMAES has good convergence capabilities. Although, some algorithms have less execution time as compared to CE-CMAES, but they have high value of R.I, mean fitness and overall costs. Overall, in terms of all the performance parameters, CE-CMAES gives the best solutions as compared to all tested algorithms. In Table 3, relative to the baseline (a last row with No LEM), all groups of agents enhance their costs/profits. In the baseline, producer agents present '0' profits due to their unwillingness to sell their energy capacity to the grid at feed-in tariffs. The reason is the CHP generators have a marginal cost higher than the feed-in tariff. As a result, generators would not sell the energy to the grid if they could not recover at least their marginal costs.
To examine the profits made by independent agents, Figure 6 shows the profits made by individual agents. Figure 6(a) indicates the individual profit/cost in the baseline scenario, which confirms that, when no LEM is available, the safest choice for generators agents is not selling anything to the grid. That fact results in zero profits/costs for them, which is not suitable for the system as a whole. On the contrary, Figure 6(b) demonstrates that the best solution found by CE-CMAES improves the condition of agents, i.e., consumer agents 1 to 3 reduce their costs; prosumer agents 4 and 6 raise their profits and agent 5 reduces its cost; and producer agents 7 to 9 increase their profits.  Finally, the convergence behaviour of the top four algorithms in terms of their fitness values have been analyzed. Figure 7 illustrates the fitness of the top four algorithms over 50,000 function evaluations. In the Figure 7, it can be seen that CE-CMAES gives slow convergence at start, but in final iterations convergence becomes fast and it achieves a   Table 4 shows, the mean fitness of CE-CMAES algorithm for different percentage usage of iterations in CE and CMAES out of total number of iterations. The comparison given in Table 4 proves that, CE-CMAES algorithm gives the best mean fitness, when both CE and CMAES are used for 50% of total iterations. This combination gives the slow convergence of CE-CMAES due to the time taken by the CMAES for the local exploitation of the search space, but finally it improves the solution in the last iterations as shown in Figure 7. Due to this reason, we used CE first for 50% of total iterations for the global exploration of search space and CMAES for remaining 50% of total iterations for the local exploitation of search space.

Percentage Contribution of CE and CMAES in CE-CMAES
However, a comparison focused purely on mean fitness, R.I., and overall cost indicates an insufficient way of comparing results. Apart from the fact that CE-CMAES outperforms all the methods in terms of above all performance parameters, it is essential to check whether all the algorithms have a statistically significant difference in fitness as given in Table 3. For that, a statistical Wilcoxon signed rank test was used [25].

Wilcoxon Signed Rank Test
The Wilcoxon signed rank test is a non-parametric test that compares two paired groups (algorithms) to detect the significant differences between their behaviors. This test measures the difference in fitness between the pairs and analyses the differences. This test was carried out in a pairwise comparison between the CE-CMAES and the rest of contestant algorithms. The Wilcoxon signed rank test is used to accept (or reject) the null hypothesis that two samples represent two separate populations [26]. In this regard, a validation of the hypothesis confirms the performance of the algorithm. In this test, the significance level is set to 5% to verify all the tested algorithms' statistical differences. The Wilcoxon signed ranks test gives a meaningful finding, if the value of `P' given by all pairwise comparison is below 0.05, and then it may be assumed that there is statistical proof that any one or more algorithms are 95% substantially different from the other in terms of fitness value. To determine whether algorithm CE-CMAES reached a statistically better solution than other tested algorithms, or if not, whether the alternative hypothesis is valid, the sizes of the ranks provided by the Wilcoxon Signed-Rank Test (i.e., T+ and T-, as defined in [25] are examined. In Table 5, '+' indicates cases in which the null hypothesis is rejected and the CE-CMAES algorithm exhibited a statistically superior performance in the pair-wise Wilcoxon Signed-Rank Test at the 95% significance level (P=0.05). As can be seen in Table  5, the Wilcoxon signed ranks test gives the value of 'P' below 0.05 for all the pairwise treatment. This implies that, CE-CMAES outperforms the competition participated algorithms in each run. In terms of R.I, mean fitness, execution time and overall cost the CE-CMAES algorithm outperforms all tested algorithms. So overall, the CE-CMAES provides the best solutions for bi-level Optimization problem.

Conclusions and Future Works
This paper has presented the application of hybrid CE-CMAES method to solve complex bi-level bidding optimization problem in the context of local energy markets. The bilevel bidding optimization is difficult to solve since all agents intent to maximize their own profits at upper level, while modifying the market-clearing price in the local market lower level as a product of their bids. Thus, it creates a strong interdependence between both levels. It was demonstrated that CE-CMAES provides the best solutions in terms of mean fitness, Ranking Index (R.I) and overall cost and profit/cost of all agents compared with the state-of-the-art algorithms that participated in the competition "Evolutionary Computation in the Energy Domain: Smart Grid Applications". In addition, Wilcoxon Signed Rank Statistical test was used to prove that CE-CMAES is statistically different from the rest of tested algorithms. For future work, the CE-CMAES could be applied to the LEM bidding problem considering the effects that such process has in the distribution system (and also considering network constraints). In addition, the LEM model could be modified to include bidding optimization in external markets as well (e.g., wholesale market, capacity markets, and ancillary services markets). Finally, new hybridized methods based in the CE-CMAES approached could be analyzed to enhance the quality of solutions.