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

Abstract

The increased penetration of renewables in power distribution networks has motivated significant interest in local energy systems. One of the main goals of local energy markets is to promote the participation of small consumers in energy transactions. Such transactions in local energy markets can be modeled as a bi-level optimization problem in which players (e.g., consumers, prosumers, or producers) at the upper level try to maximize their profits, whereas a market mechanism at the lower level maximizes the energy transacted. However, the strategic bidding in local energy markets is a complex NP-hard problem, due to its inherently nonlinear and discontinued characteristics. Thus, this article proposes the application of a hybridized Cross Entropy Covariance Matrix Adaptation Evolution Strategy (CE-CMAES) to tackle such a complex bi-level problem. The proposed CE-CMAES uses cross entropy for global exploration of search space and covariance matrix adaptation evolution strategy for local exploitation. The CE-CMAES prevents premature convergence while efficiently exploring the search space, thanks to its adaptive step-size mechanism. The performance of the algorithm is tested through simulation in a practical distribution system with renewable energy penetration. The comparative analysis shows that CE-CMAES achieves superior results concerning overall cost, mean fitness, and Ranking Index (i.e., a metric used in the competition for evaluation) compared with state-of-the-art algorithms. Wilcoxon Signed-Rank Statistical test is also applied, demonstrating that CE-CMAES results are statistically different and superior from the other tested algorithms.

1. Introduction

Smart grid technologies open the possibility of a large penetration of distributed renewable energy sources, posing several challenges for utilities and operators [1]. In this paradigm, local energy markets (LEMs) enable small agents to trade energy at the local level and contribute to the reduction of environmental pollutants [2]. Thus, LEMs empower end-users, incentivizing their organization in energy communities and the implementation of fully transactive energy systems [3]. In fact, prosumers with low production capability were unable to participate in electricity markets, due to legislative restrictions. LEMs respond to this issue by facilitating a platform for prosumers, consumers and producers, to actively participate in energy trading activities [4]. The energy transactions in LEMs can be modeled as a bi-level optimization problem in which all agents strive to maximize their profits by optimizing their bidding strategies.

The bidding optimization problem in local electricity markets is a problem that can be tackled under different perspectives and solved by different approaches. For instance, the participation of players in LEMs is typically formulated as a bi-level optimization problem and, despite taking different assumptions, modifications, and hybridizations of algorithms, finding optimal and near-optimal solutions to the problem remains a major challenge in the field.

Thus, it is important to develop a robust optimization method, which can provide near optimal solutions solving this complex problem. Therefore, this work proposes the use of a new hybrid Cross Entropy Covariance Matrix Adaptation Evolution Strategy (CE-CMAES) to address a complicated bi-level bidding optimization issue in the environment of LEMs. The main contributions of this work are as follows:

• An efficient computational intelligence method to solve an inherently non-linear and complex bidding optimization problem in LEMs.
• The problem is modeled as a multi-period bi-level optimization problem in which competitive agents at the upper level strive to maximize their profits (i.e., a multi-leader problem). The agent bids/offers modify the market clearing price response determined at the lower-level problem (i.e., a single-follower problem), resulting in a strong interdependence of their decisions.
• A detailed description of the first ranked winner optimization algorithm entitled “CE-CMAES”. The algorithm achieved the best performance in solving the testbed “bi-level optimization of end-users’ bidding strategies in local energy markets” at the international competition “Evolutionary Computation in Uncertain Environments: A Smart Grid Application” held at the Genetic and Evolutionary Computation Conference (GECCO 2020) and IEEE world Congress on Computational Intelligence (WCCI 2020).
• The effectiveness of the proposed algorithm was tested in a case study considering a power system distribution network with renewable energy sources. Also, a comparative analysis was performed showing that CE-CMAES achieves better profits for all the agents compared with different state-of-the-art algorithms.

The remaining article is organized as follows. After the introduction in Section 1, Section 2 presents a literature review of related works. The formulation of the bi-level optimization problem is discussed in Section 3. Section 4 discusses the CE-CMAES algorithm in detail. In Section 5, comparative analysis of the proposed method against state-of-the-art algorithms is carried out in a practical case study considering a distribution system. Finally, Section 6 outlines the main conclusions and future research work.

2. Literature Review

Bidding in LEMs is a problem that has been addressed considering different perspectives and utilizing a wide variety of solution methods. For instance, ref. [5] proposed a bi-level optimization problem to optimize offering prices of energy producers. The problem was later converted into a mixed integer linear programming (ILP) problem, but demand side bidding was ignored and not optimized. An optimal bidding strategy of transactive agents, dividing the formulation into two sub-problems to reduce number of variables, was proposed in [6]. Nonlinear constraints were replaced by linear inequalities to reduce the computational burden, but the effect of renewables was ignored in the problem formulation. A similar bi-level framework was used in [7] to optimize the charge and discharge power of a battery storage system at the upper level, and the Active-Reactive optimal power flow at the lower level. Similarly, the Active-Reactive reserve power flow optimization and balancing of energy curtailment in the presence of wind energy in 41 bus electricity markets was explored in [8,9].

Authors of [10] proposed a two-stage robust model to optimize the bidding in day-ahead and real time markets, considering a virtual power plant environment. The optimization problem was formulated using a linear programming (LP) model, simplifying the problem up to a level in which the quality of solutions was compromised. In [11] a hybrid optimal bidding strategy for a microgrid in a day-ahead market under the influence of intermittent distributed generation and price responsive loads was proposed. Stochastic models were used for the day-ahead market price and uncertain output of intermittent DG, while MILP was used to limit power imbalance, considering the uncertainty of the real time market price. As expected, the hourly bidding curves obtained by the linear model were less accurate than the original nonlinear problem. In [12] a robust based bidding and offering strategy was proposed for a price-maker energy storage facility participating in a day-ahead market. In [13], a unique dual bidding technique for multi-hierarchical P2P energy trading, that incorporates both intra- and inter-community trade, while accounting for uncertainties in solar irradiance and temperature, was presented. A bi-level max-min MILP approach is used to represent the participation strategy and manage the risk of uncertainty associated with forecasted hourly generation and demand price. The technique demonstrated significant reduction in costs and increase in profits.

On the other hand, some authors have explored the use of metaheuristics to solve bi-level problems in LEMs. An example of such applications can be found in [14], where Particle Swarm Optimization (PSO) is used to solve a bi-level problem in which biddings of power producers are optimized at upper level considering monthly contracts and balancing markets. The work used a standard PSO, showing poor local exploitation search ability and, thus, suffering from premature convergence. The work in [15] also used a standard Genetic Algorithm (GA) to optimize the bidding price of retailers, considering probability-based estimation of market prices. The results showed that GA is quite computationally expensive, affecting the application capabilities of the approach. In [16], an Ant Colony Optimization (ACO) algorithm was proposed to optimize the biddings of all market participants in an LEM environment. In [17] a probabilistic approach to find the optimal bidding curve to be submitted by an aggregator in day-ahead markets was proposed. Readers may refer also to [17], finding more details about two-tier LEM mechanisms. Finally, optimal day-ahead bidding strategies were analyzed in [18] from the perspective of multi-agent deep deterministic policy gradient to approximate the Nash equilibrium.

In a particular line of research, similar to the proposed work, the bi-level optimization model for bidding in LEM is solved using the vortex search single-solution-based metaheuristic algorithm, which takes into account the uncertainties of renewables and demand, in [19]. The solutions obtained by state-of-the-art evolutionary computation methods solving the bi-level bidding optimization problem in LEM are statistically compared in [20]. Later on, ref. [21] proposed the use of ACO, HyDE-DF, vortex search (VS), and estimate distribution algorithm (EDA) to solve a small variation of the problem.

The literature review reveals that the bidding optimization problem of market participants is typically formulated as a bi-level optimization problem and, despite having different assumptions, modifications, and hybridizations of algorithms, finding optimal and near-optimal solutions to the problem remains a major challenge in the field. Concerning the application of metaheuristics, the no free lunch theorem [22] of optimization also proves that designing an algorithm that consistently provides effective and robust solutions for all types of complex non-linear problems is not possible. Recently, hybridization of evolutionary algorithms has attracted attention for solving energy market problems, due to their ability to handle non-linearity and uncertainty [20,23,24]. Typically, a bi-level optimization problem is inherently highly nonlinear [25], non-convex, non-differentiable, discontinued and NP-hard [26]. Also, it is common to have a nonlinear lower-level problem with multiple optimal solutions, as shown in Figure 1. In addition, upper- and lower-level problems are mutually dependent.

The following sections introduce the optimization model and the developed evolutionary algorithm CE-CMAES, to later evaluate its performance under a realistic case study.

3. Optimization Problem Formulation

A day-ahead LEM bidding optimization problem is proposed in [27], in which agents submit offers and bids to maximize their profits and minimize their costs, respectively. Three types of agents are considered in the LEM, 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. We assume that all agents can trade energy in the LEM with prices between the feed-in tariff $\left(c^F\right)$ and the grid electricity tariff $\left(c^G\right)$. Also, it is assumed that $c^F<c^G$, and, therefore, direct energy transactions with the main grid, are less beneficial than between agents in the LEM. The transmission and distribution levels are not considered in this formulation under the assumption that agents trade energy between limits imposed by a DSO [28]. Therefore, voltage violation problems have not been considered in this work. We acknowledge that voltage violation problems for low-voltage networks can play an important role in future energy markets. Therefore, it will be considered in future work. Figure 2 shows the LEM scenario.

3.1. Bi-Level Optimization Problem Formulation

The LEM bidding optimization problem is modeled as a bi-level optimization problem in which the upper-level problem corresponds to the maximization of agents’ profits and the lower-level problem to the maximization of energy transacted in the LEM [27]. Thus, the solution of the lower-level affects the upper-level by modifying the profits of all agents.

3.1.1. Upper Level Problem

Assume a set of producer agents $J=\{1,2,\dots ,N_p\}$ and a set of consumer agents $I=\{1,2,\dots ,N_c\}$, where each agent $j$ wants to maximize its profits and each agent $i$ wants to minimize its costs.

The optimization problem (profit maximization) for the producers can be formulated as Equation (1):

$$max.P_j={\displaystyle \sum _i{c}_p{x}_{j,i}}+c^{F}Esel{l}_j^{\mathrm{grid}}-c^{mc}{G}_j$$

(1)

where $P_j$ is the profit of agent $j$, $c_p$ is the LEM clearing price (equal for buyers and sellers), $x_{j,i}$ is the energy sold by agent $j$ to agent $i$, $c^F$ is the feed-in tariff, $Esel{l}_j^{\mathrm{grid}}$ is the energy sold by agent $j$ to grid, and $c^{mc}{G}_j$ is the marginal cost associated to $j$. $c^{mc}$ is taken as zero for PV generation. Let, $c^{mc}=c_{CHP}^{mc}\left(G_j\right)$ be the marginal cost associated to a Combined Heat and Power (CHP) generator, defined as a monotonically decreasing function [29], as given in Equation (2):

$$c_{CHP}^{mc}\left(G_j\right)=\frac{b_{CHP}\sqrt{G_j}}{G_j}$$

(2)

where $b_{CHP}$ is a constant cost factor of the CHP generation unit and $G_j$ is the energy produced by the CHP unit.

On the other hand, the optimization problem (cost minimization) for consumers can be formulated as Equation (3):

$$min.C_i={\displaystyle \sum _j{c}_p{x}_{j,i}}+c^{G}Ebu{y}_j^{\mathrm{grid}}$$

(3)

where $C_i$ is the cost of agent $i$, $x_{j,i}$ is the energy sold by agent $j$ to agent $i$ in the LEM, $c^G$ is the Grid price, and $Ebu{y}_j^{\mathrm{grid}}$ is the energy bought by agent $i$ from the grid.

3.1.2. Lower Level Problem

The agents’ profits and costs depend on the LEM clearing price, which is determined in the lower-level problem and depends upon the bidding process. The lower-level problem is modeled as an asymmetric pool market in which the offers and bids are allocated using a merit order mechanism to determine the energy supply and demand curves [10]. In other words, the lower-level problem determines the clearing price $\left(cp\right)$ and energy transacted between the agents. The reader can be referred to [27] for a detailed modeling of the lower-level problem (market clearing price procedure). After that, a simple deterministic mechanism is used to determine the amount of energy traded with the grid (i.e., variables $Esel{l}_j^{\mathrm{grid}}$ and $Ebu{y}_j^{\mathrm{grid}}$ in Equation (1) and Equation (3) respectively.) Producer $j$ determines the energy sold to the grid, considering its marginal cost $c^m$, as in Equation (4):

$$Esel{l}_j^{\mathrm{grid}}=\left\{\begin{array}{ll}{P}_{max,j}-{\displaystyle \sum _i}{x}_{j,i}& if{P}_j^{\mathrm{marginal}}>P_j\\ O& Otherwise\end{array}\right\}$$

(4)

where, $P_j^{\mathrm{marginal}}$ are the profits that agent $j$ could make considering the marginal cost of producing its maximum capacity and selling the LEM non-traded energy to the grid and $P_j$ 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. For PV generation, where the marginal production cost can be considered zero, the agents will always sell their remaining energy to the grid. Consumer $i$ decides the energy bought from the grid following Equation (5):

$$Ebu{y}_j^{\mathrm{grid}}=L_{max,i}-{\displaystyle \sum _j{x}_{j,i}}$$

(5)

where, $Ebu{y}_j^{\mathrm{grid}}$ is the energy bought by consumer $i$ from the grid, $L_{max,i}$ is the maximum demand of consumer $i$, $x_{j,i}$ is the energy bought by consumer $i$ from producer $j$ in LEM. So, using Equations (4) and (5), the second terms of Equations (1) and (3) can be found.

3.1.3. Formulation of Objective Function

In Equations (1) and (3), the clearing price is the equilibrium price at which supply and demand curves intersect. The agents’ profits/costs are influenced by the LEM clearing price $c_p$, which is determined in the lower-level problem and depends on the bidding process. 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 the 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):

$$\max .{\displaystyle \sum _j\left(\frac{P_j}{N_p}\right)}-{\displaystyle \sum _i\left(\frac{C_i}{N_c}\right)}$$

(6)

where, $\left(P_j\right)$ is the profit earned by the producer agent j by selling energy to consumer agent $i$ and to the grid; $\left(C_i\right)$ is the cost of consumer agent $i$ by buying the energy from the agent j and grid. Thus, Equation (6) reflects the maximization of the overall mean profits. The profits $\left(P_j\right)$ and costs $\left(C_i\right)$ are conflicting objectives, since all agents try to achieve the best results for themselves.

3.2. Representation of Solutions and Fitness Function

As showed in Equation (6), the bi-level optimization model tries to maximize profit of all agents by executing optimal bidding of agents in the LEM. Let, $K=\{1,2,\dots ..N_k\}$ be a set including all agents (producers and consumers). The aim is to find the best tuple $\left(q_{k,}{p}_k\right)\forall k\in K$ representing the quantity of active power and optimal price to bid in the LEM for each agent. The bidding is done for all $t\in T$ optimization periods. $T=24$ periods are considered for the day-ahead market. Thus, the vector of optimization variables $\overrightarrow{x}$ is represented as $\overrightarrow{x}=\left(\left[q_{k,t}\right]{\displaystyle \cup \left[p_{k,t}\right]}\right)$, where $q_{k,t}$ is the quantity and $p_{k,t}$ is the price that the $k^{th}$ agent submits to the LEM at time period $t$. A sign convention is used to define consumer and producer types. 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 $\left[0,L_{\max }\right]$ (i.e., between 0 and their maximum consumption limit), while producer agents can submit their offers within the bounds $\left[-P{C}_{\max ,}0\right]$ (i.e., between 0 and their maximum power production capacity). The minimum and maximum boundary values for prices are in the range $\left[c^F,c^G\right]$ and are the same for all agents. Figure 3 shows the structure of a solution (i.e., a particle in our algorithm). Such particles are evaluated in the fitness function returning the mean average profit of all agents, plus the standard deviation, as per Equation (7):

$$Fitness\left(\overrightarrow{x}\right)=-mean\left(Profits\right)+std(Profits)$$

(7)

where $mean\left(profits\right)$ and $std\left(profits\right)$ 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 Equation (6) is the objective function to maximize the overall profits, where Equation (7) is the fitness function, which is the sum of the objective function and standard deviation. 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.

4. CE-CMAES Method

CE-CMAES is a sequential hybridization of Cross Entropy (CE) proposed by Rubinstein [30] and CMAES was proposed by Nikolaus Hansen [31]. It has been applied to solve the dynamic optimal power flow problem in [32]. Both are population-based iterative metaheuristic optimization methods, like Particle Swarm Optimization [33] or Differential Evolution [34]. However, unlike conventional metaheuristics, which directly act upon prospective solutions (particles) to obtain optimal ones, these methods generate “distribution” of prospective solutions updating the parameters of the distribution to obtain the optimal solutions. The CE method was used for global exploration of the search space in the initial iterative process (50% of the iterations). CE has a small number of parameters to be tuned, making it easy to be implemented. The CMAES method was used for local exploitation of search space for the remaining 50% iterations. The update of mean values in CMAES maximizes the likelihood of successful candidates, whereas the update of the Covariance matrix increases the likelihood of successful steps towards optimal solutions. Furthermore, 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 behind the CE method.

4.1. 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) $f\left(.;\theta \right)$, where $\theta$ is the vector of parameters to be optimized. Generally, $f\left(.;\theta \right)$ is a Gaussian distribution parameterized by its mean $m$ and variance ${\sigma }^2$, i.e., $\theta =\left(m,{\sigma }^2\right)$. Secondly, a threshold value $\left(\delta \right)$ of the fitness function is selected and only those particles, whose fitness values are less than the threshold value, are considered, i.e.,$f\left(x\right)<\delta$. Such particles are known as “elite particles” $\left(\mu \right)$. Then, the new parameterized distribution function $f\left(.;{\theta }^n\right)$ with elite particles is updated to coincide with the target distribution function $f\left(.;{\theta }^{\ast }\right)$ by minimizing Kullback-Leibler divergence. This procedure completes one iteration. In the subsequent generations, a family of distribution functions $f\left(.;{\theta }^{\left(1\right)}\right)\dots f\left(.;{\theta }^{*}\right)$ are produced according with ${\delta }^{\left(1\right)}\dots {\delta }^{(\ast )}$ to closely reach the optimal distribution function $f\left(.;{\theta }^{*}\right)$. The following section shows the step-by-step procedure of CE-CMAES.

Step: 1 Initialization of mean and standard deviation

Take generation $(g)=0$, total no. of generations $g_{max}=500$, number of particles $\left(N{P}_a\right)$, dimension of the problem $\left(D\right)$, elite particles $\left(\mu \right)$$\left(20\%-50\%ofN{P}_a\right)$, smoothing parameters $\alpha$ and $\beta$, $m^{\left(g\right)}\in {ℝ}^D$ is the mean value of the search distribution for each dimension at generation $g,{\sigma }^{\left(g\right)}\in {ℝ}^{+}$ is standard deviation at generation $g$. $x_{min}$ and $x_{max}$ are the minimum and maximum boundary limits of $D^{th}$ dimension particle. (i.e., $\left[0,L^{max}\right]$ for consumer bid quantities, $\left[-P{C}_{max,}0\right]$ for producer offers quantities and $\left[c^F,c^G\right]$ for bid prices). The mean and standard deviation of population are initialized as Equation (8) and Equation (9), respectively: Equation (9) is used to initialize ${\sigma }^{\left(0\right)}$.

$$m^{\left(0\right)}=\frac{\left(x_{min}+x_{max}\right)}{2}$$

(8)

$${\sigma }^{\left(0\right)}=\frac{\left(x_{max}-x_{min}\right)}{4}$$

(9)

Step: 2 Generation of population of particles

Generate the population of particles from the sampling distribution $N\left(m^{\left(g\right)},{\sigma }^2{}^{\left(g\right)}\right)$ as Equation (10):

$$x_i^{(g+1)}=m^{(g)}+{\sigma }^{(g)}randn()fori=1,\dots ..,N{P}_a$$

(10)

where $N{P}_a$ is the total number of particles (Population Size), $x_i^{(g+1)}\in {ℝ}^D$ is the i^th particle obtained at generation $\left(g+1\right)$, $m^{\left(g\right)}\in {ℝ}^D$ is the mean value of the search distribution for each dimension at generation $g$, ${\sigma }^{\left(g\right)}\in {ℝ}^{+}$ is Step-size (Standard Deviation) at generation $g$, $randn()$ is a normally distributed random variable with parameters $N\left(0,1\right)$.

Step: 3 Adjustment of search distribution limit violations

The maximum and minimum boundary limits of all particles are checked and set as Equation (11):

$$\begin{array}{l}{x}_i^{\left(g+1\right)}=x_i^{\left(g+1\right)},if{x}_{min}<x_i^{\left(g+1\right)}<x_{max}\\ x_i^{\left(g+1\right)}=x_{min},if{x}_i^{\left(g+1\right)}\le x_{min}\\ x_i^{\left(g+1\right)}=x_{max},if{x}_i^{\left(g+1\right)}\ge x_{max}\end{array}$$

(11)

Step: 4 Execution of lower-level problem

The bids and offers decoded from each particle are used to determine the market clearing price $\left(cp\right)$ through a merit order mechanism. This mechanism results in the energy transacted between the agents using Equations (4) and (5).

Step: 5 Execution of Upper-Level Problem

Calculate the profit/cost of all agents using Equations (1) and (3).

Step: 6 Fitness Function Evaluations and ranking

The fitness of all particles is evaluated with Equation (7) which calculates the mean average profit of all agents and its standard deviation. Sort (rank) all fitness values in ascending order as given in Equation (12).

$$f\left(x_1\right)<f\left(x_2\right)<\dots <f(N{P}_a)$$

(12)

where,$f(x_i)$ is the fitness of $i^{\mathrm{th}}$ particle, $x_1$ is a Global Best $\left(G_{Best}\right)$ particle having the minimum fitness among all particles. We also consider the top best 20% particles as “elite” particles.

Step: 7 Updating the mean and standard deviation of elite particles

The mean $\left(m_{\mu }^{\left(g+1\right)}\right)$ of the selected elite particles is found by Equation (13):

$$m_{\mu }^{\left(g+1\right)}=\frac{1}{\mu }{\displaystyle \sum _{i=1}^{\mu }{x}_{i:N{P}_a}^{\left(g+1\right)}}$$

(13)

where, $x_{i:N{P}_a}^{\left(g+1\right)}$ is the $i^{\mathrm{th}}$ best particle among whole population at $(g+1)$ iteration. The index $i:N{P}_a$ denotes the index of the $i^{\mathrm{th}}$ rank particle. The standard deviation $\left({\sigma }_{\mu }^{\left(g+1\right)}\right)$ of elite particles is found by Equation (14):

$${\sigma }_{\mu }^{\left(g+1\right)}=\sqrt{\frac{1}{\mu -1}{{\displaystyle \sum _{i=1}^{\mu }\left(x_{i:N{P}_a}^{\left(g+1\right)}-m^{\left(g+1\right)}\right)}}^2}$$

(14)

Step: 8 Smoothing of mean and standard deviation of population

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 shown in Equation (15):

$$m^{\left(g+1\right)}=\alpha m_{\mu }^{\left(g+1\right)}+\left(1-\alpha \right)m^{\left(g+1\right)}$$

(15)

where smoothing parameters $\alpha =0.9,\beta =0.1$:

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 than the standard deviation of all particles, which requires more exploration and thus more smoothing as given in Equation (16).

$${\sigma }^{\left(g+1\right)}=\beta {\sigma }_{\mu }^{\left(g+1\right)}+\left(1-\beta \right){\sigma }^{\left(g+1\right)}$$

(16)

Step: 9 Increment of generation count

Set $g:=g+1$

Go to step #2 and repeat steps 2–9 until 50% of $g_{max}$ are completed. So, the CE method is used for global exploration. Subsequently, the CMAES method is applied for local exploitation, which is discussed in the following section.

4.2. CMAES Method

Step: 10 Initialization of mean of each dimension

Global Best $\left(G_{Best}\right)$ particle obtained from the CE method is considered as the mean value of the corresponding dimension of CMAES, as shown in Equation (17). As a result, the new search distributions (individuals) will be generated near the sub-optimal solutions. So, the local exploitation of the search space becomes more effective than with the traditionally initialized mean value of each dimension.

$$m^g\leftarrow G_{best}$$

(17)

Step: 11 Sampling and generation of population of particles

The Covariance matrix $C$ is initialized to Identity matrix with size $\left(D*D\right)$. A population of particles (solutions or individuals) is generated by sampling a multivariate normal distribution as shown in Equation (18).

$$x_i^{(g+1)}=m^{(g)}+{\sigma }^{(g)}N\left(0,C^{(g)}\right)fori=1,\dots ..,N{P}_a$$

(18)

where, $N{P}_a$ is the total number of particles (Population Size), $x_i^{(g+1)}\in {ℝ}^D$ is the i^th particle obtained at generation $\left(g+1\right)$, $m^{\left(g\right)}\in {ℝ}^D$ is the mean value of the search distribution for each dimension at generation $g$, ${\sigma }^{\left(g\right)}\in {ℝ}^{+}$ is Step-size (standard deviation) at generation $g$, $C^{\left(g\right)}\in {ℝ}^{\left(D*D\right)}$ is the covariance matrix at generation $g$.

Step: 12 Adjustment of search distribution limit violations

The maximum and minimum boundary limits of all particles are checked and set according to Equation (11).

Step: 13 Execution of lower-level problem

Particles’ bids and offers are used to determine market clearing price $\left(cp\right)$ through merit order mechanism. Obtain energy transacted between the agents using Equations (4) and (5).

Step: 14 Execution of upper-level problem

Calculate the profit/cost of all agents using Equations (1) and (3).

Step: 15 Fitness function evaluations and ranking

The fitness of all particles is evaluated by using Equation (7) which calculates 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 $G_{Best}$ and top 50% elite particles $\left(\mu \right)$ which have better fitness than the remaining 50% particles.

CMAES is an iterative method in which $m^{\left(g+1\right)},C^{\left(g+1\right)}$ and ${\sigma }^{\left(g+1\right)}$ are updated in every generation to obtain better solutions.

Step:16 Updating the mean $\left(m^{\left(g+1\right)}\right)$

The updated mean of the search distribution is a weighted average of selected $\mu$ elite particles from the whole population according to Equation (19).

$$\begin{array}{c}{m}^{\left(g+1\right)}={\displaystyle \sum _{i=1}^{\mu }{w}_i}{x}_{i:N{P}_a}^{\left(g+1\right)}\\ m^{\left(g+1\right)}={\displaystyle \sum _{i=1}^{\mu }{w}_i}{x}_{i:N{P}_a}^{\left(g+1\right)}\\ w_1\ge w_2\ge \dots ..\ge w_{\mu }>0\end{array}$$

(19)

where, $x_{i:N{P}_a}^{\left(g+1\right)}$ is the $i^{\mathrm{th}}$ best particle among the whole population at $(g+1)$ generation. The index $i:N{P}_a$ denotes the index of the $i^{\mathrm{th}}$ rank particle. $w_1$ is the weight applied to $G_{Best}$ particle. As a higher weight is applied to $G_{Best}$ particle, the new individuals will be pulled towards it. The update of mean values maximizes the likelihood of successful candidates.

Step:17 Updating the covariance matrix $\left(C^{\left(g+1\right)}\right)$

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 for the next generation is produced from the current covariance matrix as per Equation (20).

$$C^{\left(g+1\right)}=\left(1-C_{\mu }{\displaystyle \sum w_i}\right)C^{\left(g\right)}+C_{\mu }{\displaystyle \sum _{i=1}^{N{P}_a}{w}_i{y}_{i:N{P}_a}^{\left(g+1\right)}{y}_{i:N{P}_a}^{{\left(g+1\right)}^T}}$$

(20)

where, $C_{\mu }$ is a learning rate for updating the covariance matrix, $0\le C_{\mu }\le 1$.

If $C_{\mu }=0$: No learning takes place, if $C_{\mu }=1$: No prior information is retained

$$\begin{array}{l}{w}_1\ge \dots \ge w_{\mu }>0\ge w_{\mu +1}\ge w_{N{P}_a}\\ {\displaystyle \sum _{i=1}^{\mu }{w}_i}=1\\ {\displaystyle \sum _{i=1}^{N{P}_a}{w}_i}\approx 0\end{array}$$

$$y_{i:N{P}_a}^{\left(g+1\right)}=\frac{\left(x_{i:N{P}_a}^{\left(g+1\right)}-m^{\left(g\right)}\right)}{{\sigma }^{\left(g\right)}}$$

(21)

Here, $y_{i:N{P}_a}^{\left(g+1\right)}$ is the part of rank $\mu$ update of covariance matrix. Then, $y_{i:N{P}_a}^{\left(g+1\right)}$ is updated by taking the sum of consecutive steps of distributed mean $m$ as per Equation (21).

The learning rate for updating the covariance matrix is decided by the following Equation (22).

$$C_{\mu }=\min \left(\frac{{\mu }_{eff}}{D^2},1-c_1\right)$$

(22)

$${\mu }_{eff}={\left({\displaystyle \sum _{i=1}^{\mu }{w}_i^2}\right)}^{-1}$$

(23)

where, ${\mu }_{eff}$ is the variance effective selection mass for the mean obtained from the Equation (23). $c_1\le 1-C_{\mu }$, is the learning rate for the rank-one update of the covariance matrix.

To enhance this step-size, the “evolution path” is also used. It is the sum of consecutive steps. To construct the evolution path $p_c^{\left(g+1\right)}$, exponential smoothing is applied to assign more weight to recent generations, according to Equation (24).

$$p_c^{\left(g+1\right)}=\left(1-c_c\right)p_c^{\left(g\right)}+\left(\sqrt{c_c\left(2-c_c\right){\mu }_{eff}}\right)\left(\frac{m^{\left(g+1\right)}-m^{\left(g\right)}}{{\sigma }^{\left(g\right)}}\right)$$

(24)

where, evolution path at generation $g$ is $p_c^{\left(g\right)}\in R^D$ and $p_c^{\left(0\right)}=0$. $c_c\le 1$ is the learning rate.

$\left(\sqrt{c_c\left(2-c_c\right){\mu }_{eff}}\right)$ is a normalization constant for $p_c$. ${\mu }_{eff=}{\left({\displaystyle \sum _{i=1}^{\mu }{w}_i^2}\right)}^{\left(-1\right)}$ is variance effective selection mass $1\le {\mu }_{eff}\le \mu$.

By using Equation (24), rank-one update of the covariance matrix $C^{\left(g\right)}$ is obtained according to Equation (25).

$$C^{\left(g+1\right)}=\left(1-c_1\right)C^{\left(g\right)}+c_1{p}_c^{\left(g+1\right)}{p}_c^{{\left(g+1\right)}^T}$$

(25)

$$c_1\approx \frac{2}{D^2}$$

(26)

where, $c_1$ is the learning rate find from the Equation (26).

The combination of Equations (23) and (25) yields covariance matrix adaption equation as given in Equation (27).

$$C^{\left(g+1\right)}=\left(1-c_1-c_{\mu }{\displaystyle \sum w_j}\right)C^{\left(g\right)}+c_1{p}_c^{\left(g+1\right)}{p}_c^{{\left(g+1\right)}^T}+c_{\mu }{\displaystyle \sum _{i=1}^{N{P}_a}{w}_i{y}_{i:N{P}_a}^{\left(g+1\right)}}{\left(y_{i:N{P}_a}^{\left(g+1\right)}\right)}^T$$

(27)

Step: 18 Updating the Adaptive Step-size $\left({\sigma }^{\left(g+1\right)}\right)$

The step-size is controlled by an evolution path, which is a sum of successive steps in each dimension and it is given in Equation (28).

$$p_{\sigma }^{\left(g+1\right)}=\left(1-c_{\sigma }\right)p_{\sigma }^{\left(g\right)}+\sqrt{c_{\sigma }\left(2-c_{\sigma }\right){\mu }_{eff}}{C}^{\left(g\right)}{}^{{}^{-1/2}}\left(\frac{m^{\left(g+1\right)}-m^{\left(g\right)}}{{\sigma }^g}\right)$$

(28)

where, $p_{\sigma }^{\left(g\right)}\in R^D$ is the conjugate evolution path at generation $g$ with $p_{\sigma }^{\left(0\right)}=0$.

The value $c_{\sigma }<1$ is the learning rate for conjugate evolution path, $\sqrt{c_{\sigma }\left(2-c_{\sigma }\right){\mu }_{eff}}$ is a normalization constant and $C^{\left(g\right)}{}^{{}^{-1/2}}$ is calculated as per Equation (29).

$$C^{\left(g\right)}{}^{{}^{-1/2}}=B^{\left(g\right)}{D}^{{\left(g\right)}^{-1}}{B}^{\left(g\right)}{}^{{}^T}$$

(29)

where, $B^{\left(g\right)}$ is an orthonormal basis of eigenvectors and the diagonal elements of the diagonal matrix $D^{\left(g\right)}$ are square roots of the corresponding positive eigenvalues. The Euclidian length $\left(\Vert p_{\sigma }^{\left(g+1\right)}\Vert \right)$ of evolution path is compared with its expected length $\left(E\Vert N\left(0,I\right)\Vert \right)$ to update step size as given in Equation (30).

$${\sigma }^{\left(g+1\right)}={\sigma }^{\left(g\right)}\exp \left(\frac{c_{\sigma }}{d_{\sigma }}\left(\frac{\Vert p_{\sigma }^{\left(g+1\right)}\Vert }{E\Vert N\left(0,I\right)\Vert }-1\right)\right)$$

(30)

where, ${\sigma }^{\left(g\right)}>0$.

The value $d_{\sigma }\approx 1$, is a damping parameter

$$\begin{array}{l}f\Vert p_{\sigma }^{\left(g+1\right)}\Vert >E\Vert N\left(0,I\right)\Vert \Rightarrow {\sigma }^{\left(g+1\right)}increases\\ if\Vert p_{\sigma }^{\left(g+1\right)}\Vert <E\Vert N\left(0,I\right)\Vert \Rightarrow {\sigma }^{\left(g+1\right)}decreases\\ if\Vert p_{\sigma }^{\left(g+1\right)}\Vert =E\Vert N\left(0,I\right)\Vert \Rightarrow {\sigma }^{\left(g+1\right)}remainsunchanged\end{array}$$

The step-size (standard deviation) is adaptively changed, which ultimately prevents premature convergence of CMAES [35] and improves local exploitation of the solution space.

Step: 19 Increment of generation count and Convergence check

Set $g:=g+1$

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.

5. 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 described in earlier sections, is presented.

We used the case study from the 2020 ERM competition [36] that considers nine agents, three of which are consumers, three prosumers (i.e., consumers with PV generation capabilities), and three 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 [37]. 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 of 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 b_chp = 0.18 EUR/kWh [29]. Feed-in and grid tariffs were set to c^F = 0.12 EUR/kWh and c^G = 0.28 EUR/kWh as in [27].

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” [36,38]. The MATLAB™ codes of the participating algorithms are available at http://www.gecad.isep.ipp.pt/ERM-competitions/2020-2 (accessed on 18 June 2022). The 11 algorithms included 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), CUMDANCauchy++: a Cellular EDA [39], 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 the 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 [36].

5.1. 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 to be tuned, i.e., α and β. Thus, a variety of experiments were carried out to evaluate their optimum values of the proposed CE-CMAES method. In these experiments, we fixed the population size, i.e., $N{P}_a=100$, and number of maximum iterations, i.e., (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. Tuning the smoothing parameters of CE-CMAES.

α	β	Mean Fitness
0.7	0.3	2.1128
0.75	0.25	2.1118
0.8	0.2	2.1112
0.85	0.15	2.1105
0.9	0.1	2.1088
0.95	0.05	2.1127
1	0	2.1184

.

5.2. 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. Fitness value of the tested algorithms over 20 runs.

Run	FITNESS VALUE (EUR)
	CE-CMAES	RDG3-DEEPSO	E-CBBO- CAUCHY- DEEPSO	EHL_PS_VNSO	CUMDANCAUCHY++: A CELLULAR EDA	HFEABC	DEEDA	DE-TLBO	GASAPSO	AJSO	HYDE-DF	PSO-GBP
1	2.113	2.131	2.131	2.149	2.149	2.153	2.153	2.203	2.181	2.298	3.094	4.396
2	2.109	2.125	2.131	2.149	2.150	2.153	2.172	2.329	2.232	2.326	3.077	4.518
3	2.105	2.137	2.138	2.149	2.150	2.174	2.169	2.383	2.376	2.201	3.055	4.282
4	2.103	2.119	2.137	2.149	2.149	2.151	2.160	2.325	2.272	2.315	3.037	4.433
5	2.106	2.125	2.137	2.149	2.151	2.149	2.151	2.264	2.328	2.341	2.862	4.276
6	2.110	2.125	2.143	2.144	2.149	2.159	2.159	2.206	2.295	2.351	2.956	4.332
7	2.110	2.107	2.138	2.149	2.149	2.152	2.156	2.149	2.158	2.266	2.994	4.279
8	2.110	2.119	2.137	2.149	2.149	2.170	2.149	2.137	2.297	2.239	2.954	4.578
9	2.105	2.131	2.133	2.149	2.152	2.149	2.165	2.453	2.305	2.239	3.016	4.379
10	2.117	2.125	2.136	2.149	2.149	2.158	2.158	2.482	2.249	2.301	2.992	4.397
11	2.098	2.177	2.137	2.149	2.150	2.159	2.159	2.171	2.313	2.308	3.038	4.377
12	2.103	2.137	2.132	2.149	2.150	2.180	2.153	2.182	2.234	2.243	2.953	4.348
13	2.108	2.125	2.135	2.149	2.153	2.157	2.165	2.277	2.250	2.318	3.022	4.362
14	2.111	2.198	2.126	2.149	2.149	2.157	2.158	2.281	2.224	2.271	3.200	4.542
15	2.103	2.184	2.136	2.149	2.149	2.171	2.164	2.212	2.307	2.255	2.952	4.473
16	2.123	2.119	2.132	2.149	2.149	2.150	2.155	2.113	2.326	2.207	3.101	4.348
17	2.119	2.131	2.125	2.149	2.149	2.167	2.149	2.119	2.206	2.260	3.140	4.139
18	2.105	2.119	2.126	2.149	2.149	2.150	2.157	2.115	2.284	2.346	3.050	4.575
19	2.101	2.137	2.138	2.149	2.150	2.172	2.172	2.176	2.336	2.327	3.047	4.238
20	2.108	2.125	2.145	2.149	2.149	2.200	2.161	2.127	2.284	2.272	3.103	4.211

. The fitness value displayed in Table 2 is the average value of the fitness function over the 50,000 function evaluations of each run. Table 2 clearly shows that CE-CMAES achieved the best fitness value in eighteen runs out of 20 runs, where RDG3-DEEPSO had the best fitness value in two out of 20 runs. The results suggest the robustness of the CE-CMAES algorithm relative to other tested algorithms. Table 2 also reveals that the eleventh run of CE-CMAES offered the best fitness value with 2.0988. It means that, in this run, CE-CMAES provided 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 are presented in Appendix A (

Table A1. Solution of CE-CMAES.

Agents		Consumer (1 to 3) Power Bid (kW)			Prosumer (4 to 6) Power Bid (kW)			Producer (7 to 9) Power Bid (kW)			Consumer (1 to 3) Price Bid (EUR/kWh)			Prosumer (4 to 6) Price Bid (EUR/kWh)			Producer (7 to 9) Price Bid (EUR/kWh)
Hours	1	0.20	0.87	0.16	0.23	0.78	0.19	−2.00	−1.26	−2.00	0.28	0.25	0.28	0.28	0.28	0.28	0.16	0.12	0.12
	2	0.12	2.07	0.08	0.10	2.88	0.09	−2.00	−2.00	−2.00	0.28	0.28	0.28	0.28	0.28	0.28	0.26	0.28	0.12
	3	0.53	1.76	0.08	0.52	1.93	0.10	−2.00	−2.00	−2.00	0.28	0.28	0.28	0.28	0.28	0.28	0.12	0.18	0.23
	4	0.11	0.08	0.08	0.09	0.06	0.09	−2.00	−0.35	−0.81	0.28	0.28	0.28	0.28	0.28	0.28	0.12	0.20	0.12
	5	0.04	0.06	0.09	0.04	0.05	0.09	−1.69	−1.10	−1.42	0.28	0.28	0.28	0.28	0.28	0.28	0.12	0.12	0.28
	6	0.05	0.08	0.08	0.04	0.10	0.07	−1.58	−1.68	−1.57	0.19	0.28	0.28	0.28	0.28	0.28	0.12	0.12	0.28
	7	0.08	0.07	0.07	0.06	0.05	0.07	−1.10	−1.84	−0.93	0.28	0.28	0.28	0.28	0.28	0.28	0.12	0.12	0.12
	8	0.00	0.97	0.13	−0.06	0.77	0.09	−1.35	−2.00	−0.84	0.23	0.28	0.28	0.28	0.28	0.28	0.27	0.12	0.28
	9	0.13	2.21	0.18	0.00	2.09	0.00	−2.00	−2.00	−2.00	0.28	0.28	0.28	0.12	0.28	0.12	0.12	0.26	0.28
	10	0.12	1.34	0.20	−0.43	0.92	−0.27	−1.34	−1.93	−2.00	0.28	0.28	0.28	0.12	0.28	0.12	0.12	0.12	0.12
	11	0.18	0.06	0.11	−0.65	−0.09	−0.77	−1.64	−0.56	−1.42	0.28	0.28	0.28	0.28	0.12	0.12	0.28	0.12	0.28
	12	0.13	0.07	0.07	−0.16	−0.55	−0.79	−0.33	−1.39	−0.56	0.28	0.28	0.28	0.19	0.27	0.12	0.28	0.12	0.12
	13	0.18	0.07	0.11	−0.67	−0.70	−0.73	−1.39	−1.05	0.00	0.28	0.28	0.28	0.28	0.12	0.28	0.28	0.28	0.28
	14	0.16	0.07	0.19	−0.92	−0.74	−0.23	−1.66	0.00	−1.13	0.28	0.28	0.28	0.12	0.12	0.18	0.28	0.12	0.28
	15	0.19	0.10	0.19	−0.40	−0.61	−0.51	−1.08	−0.98	−0.80	0.28	0.28	0.28	0.28	0.12	0.12	0.28	0.12	0.12
	16	0.16	0.07	0.29	−0.51	−0.52	−0.16	−0.94	−0.79	−1.00	0.28	0.28	0.28	0.12	0.12	0.12	0.15	0.28	0.12
	17	0.17	0.05	0.31	0.00	−0.30	0.00	−2.00	0.00	−0.17	0.28	0.14	0.28	0.13	0.28	0.28	0.12	0.12	0.12
	18	0.12	0.18	0.23	−0.28	0.00	0.00	−1.87	−2.00	−0.78	0.28	0.28	0.28	0.28	0.28	0.28	0.12	0.21	0.25
	19	0.16	0.19	0.25	0.00	0.00	0.02	−2.00	0.00	−0.99	0.28	0.28	0.28	0.21	0.18	0.20	0.28	0.12	0.21
	20	0.29	0.67	0.22	0.17	0.44	0.12	−2.00	−1.04	−1.09	0.28	0.28	0.28	0.28	0.28	0.28	0.28	0.27	0.12
	21	0.42	0.38	0.26	0.38	0.30	0.25	−2.00	−1.02	−1.29	0.28	0.28	0.28	0.28	0.28	0.28	0.23	0.12	0.12
	22	0.32	0.39	0.21	0.33	0.28	0.24	−2.00	−0.89	−0.83	0.28	0.28	0.28	0.28	0.22	0.28	0.12	0.14	0.17
	23	0.22	1.55	0.24	0.19	1.45	0.22	−2.00	−2.00	−1.72	0.28	0.28	0.28	0.28	0.28	0.28	0.15	0.28	0.14
	24	0.26	0.24	0.26	0.31	0.26	0.29	−1.05	−0.71	−2.00	0.28	0.28	0.28	0.28	0.26	0.28	0.12	0.28	0.28

).

All algorithms were run 20 times, and the results were recorded to obtain the mean value of those 20 runs.

Table 3. Overall costs/profits, mean fitness, standard deviation, execution time and R.I achieved by the tested algorithms.

Algorithms	Overall Costs (EUR)	Costs/Profits by Group of Agents (EUR)			MeanFit (EUR)	Std	R.I	Average Time for 20 Runs (min)
		Costs		Profits
		Consumers	Prosumers	Producers
CE-CMAES	3.1840	2.7850	0.2950	0.1040	2.1088	0.0061	2.1148	17.0210
RDG3-DEEPSO	3.2528	2.8038	0.2978	0.1513	2.1353	0.0236	2.1589	16.5985
E- CBBO_Cauchy_DEEPSO	3.2762	2.8286	0.3191	0.1285	2.1536	0.0402	2.1938	26.2037
EHL_PS_VNSO	3.2732	2.7954	0.2889	0.1888	2.1446	0.0045	2.1491	25.4820
CUMDANCauchy++: a Cellular EDA	3.2867	2.7984	0.2896	0.1987	2.1500	0.0006	2.1506	13.8346
HFEABC	3.3081	2.8066	0.3018	0.1997	2.1584	0.0058	2.1642	1.8408
DEEDA	3.3125	2.8152	0.2974	0.1999	2.1605	0.0092	2.1697	16.9520
DE-TLBO	3.4716	2.9555	0.3957	0.1205	2.2357	0.1115	2.3473	16.0524
GASAPSO	3.5583	3.0624	0.4131	0.0828	2.2734	0.0554	2.3288	1.8695
AJSO	3.6664	3.0856	0.4606	0.1202	2.3130	0.0797	2.3927	3.6674
HyDE-DF	4.9823	3.9897	1.1648	−0.1722	3.0326	0.0773	3.1099	1.8578
PSO-GBP	7.8080	5.3994	2.1901	0.2185	4.3410	0.1493	4.4903	16.8510
No LEM (baseline)	8.9985	6.1541	2.8444	0	4.9705	0	4.9705	0

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 average 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 could attain in the absence of an 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. Also notice that, CE-CMAES achieved the best overall costs of 3.1840 EUR for the system compared to all tested algorithms. In terms of execution time, CE-CMAES took an average of 17.0210 min to solve the problem, which was lower than the EHL_PS_VNSO (25.4820 min). It proves the advantage of hybridization of CMAES with CE, as the best solution obtained by CE was used for the initialization of CMAES. As a result, CMAES was near to the sub-optimal solutions. Therefore, CE-CMAES has good convergence capabilities. Despite some algorithms having faster execution times as compared to CE-CMAES, they had worse values of R.I, mean fitness and overall costs. Overall, the CE-CMAES results suggest it provided the best solution to this testbed, as compared to all tested algorithms.

In Table 3, compared with the baseline (the last row “No LEM”), all groups of agents enhanced their costs/profits. In the baseline, producer agents presented ‘0’ profits, due to their unwillingness to sell their energy capacity to the grid at the feed-in tariff. The reason was that CHP generators had a marginal cost higher than the feed-in tariff and were not eligible for the incentive. Even if eligible, these 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 the involved agents, Figure 6 shows the profit/cost made by each agent. In Figure 6a,b, a positive value indicated cost, and a negative value indicated the agents’ profits. Here, the profit and cost of prosumer/small generators were based on the difference between marginal cost and feed-in/grid tariffs. Figure 6a shows the individual profit/cost in the baseline scenario, which confirmed that, when LEM was unavailable, the safest choice for generators’ agents was not selling to the grid (see bars 7, 8, and 9). That situation resulted in zero profits/costs for producers, which was not suitable without bilateral contracts. Also, since there was no LEM in the base scenario, prosumers 4 and 5 were unable to sell their excess of PV energy and this increased the cost of prosumer 6. Figure 6b demonstrates that the best solution found by CE-CMAES improved the condition of all agents in terms of costs and profits when the LEM was available. As a result, consumer agents 1 to 3 reduced their costs; prosumer agents 4 and 6 increased their profits, and prosumer 5 reduced its cost; and finally producer agents 7 to 9 increased their profits.

Finally, the convergence behavior of the top four algorithms, in terms of their fitness values, was analyzed. Figure 7 illustrates the fitness of the top four algorithms over 50,000 function evaluations. In Figure 7, CE-CMAES gave slow convergence at the start. Nevertheless, convergence became fast in the final iterations, and it achieved a better final solution than all evaluated algorithms. The CE-CMAES had a slow convergence due to the 50–50% contribution of CE and CMAES algorithms.

5.3. Percentage Contribution of CE and CMAES in CE-CMAES

Table 4. Mean Fitness of CE-CMAES algorithm for a different contribution of iterations (%) in CE and CMAES.

CE	CMAES	Mean Fitness
10%	90%	2.8096
20%	80%	2.7466
30%	70%	2.7546
40%	60%	2.1107
50%	50%	2.1088
60%	40%	2.1132
70%	30%	2.119
80%	20%	2.1823
90%	10%	2.3354

depicts the mean fitness of the CE-CMAES algorithm for a different contribution of iterations in CE and CMAES out of the total number of maximum iterations, expressed in percentage. The comparison given in Table 4 proves that CE-CMAES algorithm achieved the best mean fitness when both CE and CMAES were balanced with 50% each of the total maximum number of iterations allowed. This combination gave the slowest convergence of CE-CMAES, due to the time taken by the CMAES for local exploitation of the search space. CMAES finally improved the solution in the last iterations, as shown in Figure 7. For this reason, we used CE first for 50% of total iterations for the global exploration of search space and CMAES for the remaining 50% of total iterations for the local exploitation of search space.

5.4. Wilcoxon Signed Rank Test

A comparison focused on mean fitness, R.I., and overall cost is insufficient to conclude the performance of the results. Apart from the fact that CE-CMAES outperformed all the methods in terms of the provided performance parameters, it was essential to check whether all the algorithms had a statistically significant difference, as given in Table 3. A statistical Wilcoxon signed rank test was used for this purpose [40].

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 the contestant algorithms. The Wilcoxon signed rank test is used to accept (or reject) the null hypothesis that two samples representing two separate populations [41]. In this regard, a validation of the hypothesis confirmed the performance of the algorithm. The significance level was set to 5% to verify all the tested algorithms’ statistical differences. The Wilcoxon signed ranks test provided a meaningful finding if the value of ‘P’ given by the pairwise comparison was below 0.05, and then it might be assumed that there was statistical evidence that the algorithm was significantly better 95% of the time.

To determine whether algorithm CE-CMAES reached a statistically better solution than other tested algorithms, or, if not, whether the alternative hypothesis was valid, the sizes of the ranks provided by the Wilcoxon Signed-Rank Test (i.e., T+ and T−, as defined in [40] were examined. In

Table 5. Wilcoxon Signed Ranks Test Results.

PAIR	p-Value	T+	T−	Winner
CE-CMAES-E-CBBO_Cauchy_DEEPSO	0	0	210	+
CE-CMAES-EHL_PS_VNSO	0	0	210	+
CE-CMAES-CUMDANCauchy++: a Cellular EDA	0	0	210	+
CE-CMAES-HFEABC	0	0	210	+
CE-CMAES-DEEDA	0	0	210	+
CE-CMAES-DE-TLBO	0.006	0	210	+
CE-CMAES-GASAPSO	0	0	210	+
CE-CMAES-AJSO	0	0	210	+
CE-CMAES-HyDE-DF	0	0	210	+
CE-CMAES-PSO-GBP	0	0	210	+

, ‘+’ indicates cases in which the null hypothesis was rejected, and the CE-CMAES algorithm exhibited statistically superior performance in the pair-wise Wilcoxon Signed-Rank Test at the 95% significance level (p = 0.05). As shown in Table 5, the Wilcoxon signed ranks test gave the value of ‘P’ below 0.05 for all the pairwise comparisons. This implied that CE-CMAES outperformed the contestant algorithms in each run. In terms of R.I, mean fitness, execution time, and overall cost, the CE-CMAES algorithm outperformed all tested algorithms. So, overall, the CE-CMAES provided the best solution for bi-level optimization problem with statistical proof.

6. Conclusions and Future Works

This paper presented the application of the hybrid CE-CMAES method to solve the complex bi-level bidding optimization problem in local energy markets. The bi-level bidding optimization is difficult to solve since all agents intend to maximize their profit at the 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 solution 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 of the competition “Evolutionary Computation in the Energy Domain: Smart Grid Applications”. Using the proposed optimization algorithm, the profit of all the producer agents was maximized by optimal scheduling of generation in LEM at optimal bidding price. It also reduced the cost of the consumer agents by buying the energy from the LEM and grid at competitive price. To confirm the superiority, Wilcoxon Signed Rank Statistical test was used to prove that CE-CMAES was statistically different from the rest of the contestants.

Regarding possible implementation of an LEM, in this work it was assumed that the various agents (consumers, producers and prosumers) and market operator posed adequate infrastructure in terms of communication and information technologies. Thus, from the obtained optimized results, the local market operator (e.g., an aggregator) might send the command signals to various market participants to change their output power to maximize their profits in a smart grid system. Our work was proposed under the idea that agents exchange energy between the limits imposed by a DSO. As a result, voltage violation issues were not addressed in this work. While voltage violation issues in low-voltage networks might be a critical aspect to be tackled in future energy markets, we recognize that it is a limitation of this work to be taken into account in upcoming work.

For future work, the CE-CMAES could be applied to the LEM bidding problem considering the distribution system constraints and also the cost of the network in the LEM. The proposed method could be additionally applied in reactive power bidding, considering active-reactive power optimization. The network states and interconnection with upper levels of the energy chain (like the distribution level) are important research avenues for future research. 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 on the CE-CMAES approached, could be analyzed to enhance the quality of solutions.