Modified Social Group Optimization to Solve the Problem of Economic Emission Dispatch with the Incorporation of Wind Power

: Economic dispatch, emission dispatch, or their combination (EcD, EmD, EED) are essential issues in power systems optimization that focus on optimizing the efficient and sustainable use of energy resources to meet power demand. A new algorithm is proposed in this article to solve the dispatch problems with/without considering wind units. It is based on the Social Group Optimization (SGO) algorithm, but some features related to the selection and update of heuristics used to generate new solutions are changed. By applying the highly disruptive polynomial operator (HDP) and by generating sequences of random and chaotic numbers, the perturbation of the vectors composing the heuristics is achieved in our Modified Social Group Optimization (MSGO). Its effectiveness was investigated in 10-unit and 40-unit power systems, considering valve-point effects, transmission line losses, and inclusion of wind-based sources, implemented in four case studies. The results obtained for the 10-unit system indicate a very good MSGO performance, in terms of cost and emissions. The average cost reduction of MSGO compared to SGO is 368.1 $/h, 416.7 $/h, and 525.0 $/h for the 40-unit systems. The inclusion of wind units leads to 10% reduction in cost and 45% in emissions. Our modifications to MSGO lead to better convergence and higher-quality solutions than SGO or other competing algorithms.


Introduction
The power generation sector has been undergoing continuous development in recent years, with a focus on diversification of energy sources and production technologies, but also on efficiency and constant innovation.Decision makers aim through regulations and policies to create a competitive and attractive framework for investors while ensuring the sustainable development of the energy system [1].In order for electricity producers and companies to remain in a competitive market, it is necessary for them to operate with the lowest possible costs and to use sources/technologies with the lowest possible environmental impact.One way to optimize operating costs while taking emissions into account is to solve the economic emission dispatch (EED) problem [2].The goal of the EED problem is to determine the optimal operating mode of energy sources to minimize the two objectives-cost and emissions-considering a given power demand, as well as the operating limits of the generating units [3].If the EED problem aims only at cost optimization (without considering emissions), then this is called the economic dispatch problem (EcD) [4].If EED only aims to optimize emissions (without considering costs), then it is called the emission dispatch (EmD) problem [5].For a more complete approach, both the economic aspects and the level of emissions released into the atmosphere must be taken into account, which implies solving the EED problem.
Initially, the solution to EcD, EmD, and EED problems focused mainly on conventional unit power plants (firing coal, gas, etc.), as they had a very large share of the power system structure.However, last year's warnings about the environmental impact of these types of power plants led to an increased use of both high-efficiency co-generation units (CU) (that simultaneously produce both electricity and heat), and renewable energy power generation systems (RESPSs) to steer the electricity generation sector towards sustainable development [6,7].The EED problem that includes CU aims at the optimal scheduling (in terms of cost and emissions) of power-only units, heat-only units, and their combination: co-generation units [8].It is a similar situation for EcD [9] or EmD [8] problems, with both aiming to manage all categories of units operating in the power system: power-only, heatonly, and CU units.For systems that include CU units, in addition to costs and constraints specific to power-only units, the EED optimization model must also consider the fuel costs related to heat-only and CU units, power-heat dependencies, and the produced-demand heat balance [10].In the case of including RESPSs in the power system, an important option is the use of wind energy, which can be converted into electricity without producing greenhouse gases.Wind energy is considered a renewable and clean energy source, having only a low secondary impact on the environment due to the manufacturing process of the equipment and its transport.Increasing the share of wind energy in the energy mix has a positive impact on the quality of the environment and reduces dependence on conventional sources.Since wind is an uncontrollable and variable source of energy, the energy production from wind units is influenced by weather conditions and wind speed.Thus, the operation of power systems that include wind units must also take into account the unpredictable variations in the power of these units [4].
The integration of wind units into EcD, EmD, or EED problems brings challenges related to wind speed distribution, mathematical optimization models, and solution algorithms.One of the most common parametric distributions used in wind speed modeling is the Weibull distribution [11][12][13], but other distributions can be considered, such as [14] Gamma, inverse Gamma, Gaussian, Burr, Halphen, etc.
In the following, some articles are presented which propose mathematical models and solution algorithms for solving EcD, EmD, or EED problems considering the uncertainties caused by unpredictable wind speed fluctuations.Thus, in [4], an optimization model is presented that includes the costs related to the overestimation and underestimation of the wind power.The case study considers a Weibull distribution of the wind speed, and it shows that the optimal solutions to the EcD problem can be influenced by factors associated with the overestimation and underestimation of wind power.Based on a linear relationship between wind speed and wind power, and considering a Weibull distribution for wind speed, in [11], an optimization model is developed for the minimization of emissions, in which the uncertainty of the wind power is included in the constraints of the model.In [12], both the cost of emissions and the cost of overestimation and underestimation of wind power are included in the bi-objective EED problem.Starting from the classical formulation of EcD and EmD problems, in [2], the objective functions related to costs and emissions are extended with terms that include the uncertainty of the wind power.The resulting models are tested on a system consisting of two conventional units and two wind units, considering a mixed Gamma-Weibull distribution for the wind speed.To solve the EED problem with the inclusion of wind power, in [15], the Honey Bee Mating Optimizer (HBMO) is applied, where a so-called 2m-point model is used to estimate the uncertainty of the wind power.In [16], a new evolutionary technique called Lightning Flash Algorithm (LFA) is proposed to solve the bi-objective EED problem by considering different levels of wind power penetration and multiple fuel sources.The LFA technique is efficient for the EED problem, having a good convergence and achieving lower costs and emissions than other algorithms.A robust and efficient optimization model for dispatching windthermal power under uncertainties is developed in [1], taking into account robustness, economics, and environmental aspects.In [17], another robust model for the optimization of the EcD problem is presented.It is based on the identification of a set of discrete bad scenarios, where the objective function aims at minimizing all the penalties associated with the bad scenarios.
The EcD, EmD, or EED problems that include or do not include wind power uncertainty are non-convex with multiple local optima and nonlinear constraints [12,15,16] requiring advanced solution algorithms.
For the study of both systems (the ones comprising conventional-only units and the ones comprising thermal-wind units), several algorithms are presented in the following lines.Considering the first category of systems, in [18], an Improved Class Topper Optimization (ICTO) has been developed by including in the classical CTO three new concepts: adaptive improvement factor, adaptive acceleration coefficient, and chaos local search, which help to enhance the exploration and exploitation capability of the ICTO.The ICTO is tested on five systems with conventional units, and it performs better, in terms of cost and emissions, than the original CTO and several other competing algorithms.A quasi-oppositional-based political optimizer is used, in [19], to solve the bi-objective EED problem considering valve-point loading effect, transmission line losses, and other constraints related to conventional units.Recently, two new population-based algorithms (Criminal Search Optimization Algorithm (CSOA) [20] and Kho-Kho optimization algorithm (KKO) [21]) have been developed and applied to solve the EED problem.Both algorithms (CSOA and KKO) show good exploitation and exploration capabilities, and they can be used to solve complex problems in various domains.Another option for solving the EED problem is presented in [22], where two metaheuristic algorithms, Exchange Market Algorithm (EMA) and Adaptive Inertia Weight Particle Swarm Optimization (AI-WPSO), are combined to obtain a new algorithm with improved global and local search abilities.The constraints of the EED problem are maintained using the multiple constraints ranking technique.The best compromise solution (BCS) obtained using EMA-AIWPSO dominates the BCSs obtained using other algorithms (such as KKO, ISA, or GSA).Also, several multi-objective algorithms have been proposed for solving the EED problem, such as multi-objective SSA [23], multi-objective cultural algorithm [24], NSGA-III algorithm [25], or multi-objective quasi-oppositional TLBO [26], each of which uses a basic metaheuristic algorithm (SSA, CA, GA, or TLBO) that is endowed with the Pareto-dominance principle to generate successive Pareto fronts.The multi-objective algorithms mentioned [23][24][25][26] have been tested on various power systems with conventional units, their performance being superior to other recognized multi-objective techniques (MODE, NSGA II, or SPEA-2).
Over the time, when trying to solve the EED problem, various techniques have been used to optimally program the units of the thermal-wind systems.Thus, in [6], the techniques of weighted goal programming and the progressive bounded constraint method are combined to generate a set of Pareto solutions that are efficient in the cost-emission space, and then extract the best compromise solution.For instance, four cases are analyzed to quantify and demonstrate the benefits of including wind units in the structure of a power system comprising only thermal units.In [12], the Artificial Bee Colony (ABC) algorithm is strengthened by including the best solution in the update equations, and in [27], chaos is inserted into the sine-cosine algorithm to obtain better-quality solutions to the EED problem.Additionally, in [13]-where wind units are considered-a number of eight metaheuristic algorithms (Flower Pollination Algorithm (FPA), Mine Blast Algorithm (MBA), Backtracking Search Algorithm (BSA), Symbiotic Organisms Search (SOS), Ant Lion Optimizer (ALO), Moth-Flame Optimization (MFO), Stochastic Fractal Search (SFS), and Lightning Search Algorithm (LSA)) are applied to identify the best solutions to the EcD, EmD, or EED problem.The FPA and BSA algorithms provide the best scheduling of thermal/wind units for cost and/or emission minimization.To study the behavior of systems in which thermal, wind, and solar units operate, in [28], the Dragonfly Algorithm (DA) is proposed for the EcD problem.Uncertainties related to the power generated by using wind and solar energy are modeled using the 2-m point estimation technique.The DA algorithm outperforms other algorithms, such as Crow Search Algorithm (CSA), Ant Lion Optimizer (ALO), oppositional RCCRO, Biogeography-Based Optimization (BBO), PSO, and GA, in terms of cost and execution time.
SGO is a metaheuristic algorithm proposed relatively recently by Satapathy SC (2016) [29], which is based on the fact that a social group of individuals has a greater ability to solve a real-life problem than a single individual.For some mathematical functions, SGO is more efficient than other well-known metaheuristic algorithms (such as [29]: DE, PSO, ABC, GA, TLBO, etc.), being an easy-to-implement algorithm, having two main phases (improving phase and acquiring phase) and a single specific parameter.However, for some mathematical functions or applications, it can provide a low performance due to an imbalance between exploration and exploitation, which ultimately leads to getting stuck in a local optimum [30].Also, SGO may suffer due to the reduced diversity of the population [31].To overcome the mentioned shortcomings and to obtain better quality solutions, there have been attempts to improve the performance of SGO by different strategies, such as modifying the relations for updating the solutions [30], inertia weight strategies [32], and hybridizations with other algorithms [31].
In this paper, some changes have been made to the original SGO algorithm, aiming to increase efficiency, resulting in the MSGO algorithm.In order to improve the exploration-exploitation balance, the commutation condition between the equations for updating the solutions in the acquiring phase was changed.Also, to increase the diversity of the population, a logistic map and a highly disruptive polynomial operator were inserted.
The main contributions of this paper are: • We propose a modified version of SGO (called MSGO) in which the way of updating and adapting the individuals in the social group is changed by inserting chaos and an HDP operator (in the original SGO only uniformly randomly generated number sequences are used).The operators associated with chaos and HDP aim at increasing the efficiency of the MSGO algorithm by reducing the number of close solutions and overcoming some drawbacks related to slow convergence.To the best of the authors' knowledge the HDP operator has never been used to improve the performance of the SGO algorithm.

•
Implementation of MSGO to solve EcD, EmD, and EED problems with or without consideration of wind units.

•
Conducting experiments to evaluate and statistically compare the effectiveness of MSGO with SGO and other well-known algorithms (or their varieties) for thermal or wind-thermal power systems of different sizes and characteristics.

Statement of the EcD Problem
The EcD problem for a power system that includes thermal and wind units aims to determine the power outputs of these generating units so that the operating cost of the entire system is minimized, and a number of technical constraints are met.For the formulation of the mathematical optimization model, we consider a power system comprising Nt thermal units and Nw wind units, and the total power demand (P D ) of the consumers in the system is considered known and constant for the period of analysis.The variables to be optimized are continuos, being represented by the output power vectors of the thermal units PT = [PT 1 , PT 2 ,. .., PT i ,. .., PT Nt ] and of wind units PW = [PW 1 , PW 2 ,. .., PW j ,. .., PW Nw ].
The objective function C(PT, PW) is represented by two components; one related to the fuel cost of the thermal units, C T (PT), and the other related to the operating cost of the wind units, C W (PW) [4]: where, C T i (PT i ) is the fuel cost corresponding to thermal unit i, and C W j (PW j ) is the operating cost corresponding to the wind unit j.
In this paper, the relationship between the fuel cost C T i (PT i ) and output power PT i of a thermal unit i is modeled via a non-convex function consisting of a quadratic and a sinusoidal term [33]: The operating cost of wind unit j, C W j (PW j ), consists of three terms [4]: the direct operating cost of unit j (c d j • PW j ), the cost corresponding to the overestimation of the wind power (c o j •E o j PW j ), and the cost corresponding the underestimation the wind power (c u j •E u j PW j ).
The wind power overestimation for unit j occurs if the estimated wind power for this unit PW j is higher than the available power represented by a random variable W j .In this situation, the difference between the two powers is covered by a reserve source and implies the reserve cost c o j •E o j PW j .In the opposite situation of the underestimation (if PW j is less than W j ), a part of the available power of unit j remains unused, which implies a penalty cost c u j •E u j PW j .Since W j is a random variable, the powers E o j PW j and E u j PW j will include the uncertainty of the wind power, and their calculation method is presented below.
For the calculation of the average powers E o j PW j and E u j PW j , it is considered that the wind speed is modeled by a Weibull distribution whose probability density function (pdf) is given by relation (5) [11]: where, V is the random variable wind speed, v denotes the wind speed (a value of V), and f V (v) is the pdf of the variable V, k is the shape parameter, and c is the scale parameter.Also, between the random variable wind power (W) and the random variable wind speed (V), we consider a linear relationship expressed as follows [11]: The use of the linear model (6) requires knowledge of three limit speeds: cut-in wind speed (v in ), rated wind speed (v r ), and cut-out wind speed (v out ). PWr is the rated output power of the wind unit, which corresponds to the speed v r .
The pdf for wind power (W) over the continuous interval (0, PWr) has the expression: where R = (v r − v in )/v in .Given relation (6) and the Weibull distribution of the wind speed, the event probabilities W = 0 and W = PWr are calculated using relations [4]: We must say that, in general, the distributions of the variables V and W, as well as the pdfs derived from them, differ depending on the location of the wind turbine units.Thus, for wind unit j, the average powers associated with the overestimation and the underestimation of the wind power are mathematically expressed as follows [4]: where f W j (w) is probability density function for the random variable W j having the form given by relation (7), while w are the values of the continuous random variable W j .The discrete probabilities Prob W j = 0 and Prob(W j = PW rj ), for unit j, are calculated using relations ( 8) and ( 9), where PWr j represents rated wind power for unit j.
In Table A1 are presented the steps to evaluate the cost related to wind power, as well as a calculation example for a wind unit.
The feasible space of EcD problem solutions is limited by the following constraints [4,12]: 1.The thermal units i must be operated between the minimum capacity (PT min,i ) and the maximum capacity (PT max,i ): If the power of the thermal unit i in the previous hour is known or specified (PT 0 i ), then the operating limits are given by the constraint: where DR i and UR i are the down-ramp and up-ramp limits of the unit i.
2. The wind units j must be operated between the minimum (PW min,j ) and maximum (PW max,j ) capacity: PW min,j ≤ PW j ≤ PW max,j , j = 1, 2,. .., Nw 3. The actual powers generated by the thermal and wind units must cover the power consumed in the system: where P D is the load demand of the system, and P L represents the transmission line losses, which can be determined considering the constant B coefficients formula: where B ij is an element of the loss coefficients matrix, B 0i is i element of the loss coefficients vector, and B 00 is the loss coefficient constant; PTW = [PTW 1 , PTW 2 ,. . ., PTW i ,. . ., PTW Nt+Nw ] is a vector that combines the powers of thermal PT and wind PW units, while PTW i represents the i th component of the PTW vector.

Statement of the EmD Problem
The EmD problem has a high level of similarity with the EcD problem, aiming to determine the PT and PW vectors, so that the emissions released into the atmosphere is minimal while maintaining some technical constrains at the units and system level.Thus, in the EmD problem, the variables to be optimized (PT and PW vectors), as well as the constraint relations ( 12)- (15), are identical to those in the EcD problem.The objective function is represented by the total amount of emissions released into the atmosphere, mathematically defined by the following relation [2]: where E T i (PT i ) and E T (PT) represent the pollutant emissions released into the atmosphere due to the operation of thermal unit i, respective of all thermal units in the analyzed power system.E W j (PW j ) and E W (PW) are the pollutant emissions produced due to the need to use other thermal units to cover the uncertainty of the availability of wind unit j, respective of all wind units considered.Average powers, E o j PW j and E u j PW j , are calculated with relations ( 10) and (11), which include the uncertainty of wind power in the estimation of pollutant emissions.
The term e o j •E o j PW j represents the emissions released into the environment due to the need to use some thermal units in the system to cover the difference between the scheduled power of the wind unit PW j (power that cannot be realized due to the unavailability of the wind resource) and its available power W j (the case of overestimation of wind power).
The term e u j •E u j PW j represents the emissions released into the environment by other thermal units due to the non-use of the full available power of wind unit j (the case of underestimation of wind power).The quantity of emissions E T i (PT i ) released into the atmosphere by a thermal unit i can be defined by relation [2]:

Statement of the EED Problem
The EED problem is similar to the EcD and EmD problems, in that the variables to be optimized and the constraints of the problem are the same, but the objective function is different.In the EED problem, the objective function Φ(PT, PW) can be formed by the weighted and normalized summation of the cost C(PT, PW) and emissions E(PT, PW) objectives [34]:  [35].Thus, the solution for which a minimum cost is obtained in the EcD problem will be considered to result in maximum emissions, and vice versa.
In current practice, system operators or decision makers look for a single solution that takes into account both objectives, which is generally the best compromise solution.By solving the EED problem for different values of the ω factor, one can estimate the set of non-dominant solutions that form the discrete Pareto front [5] (from which the BCS between the two objectives is extracted).To extract the BCS from the set of solutions belonging to the Pareto front, a fuzzy approach is used [5].Thus, the solutions in the Pareto front are ordered according to one of the objectives, then for each objective i and non-dominant solution s, a fuzzy membership function µ i,s is assigned.This is defined by relation [5]: µ i,s = (f i,max − f i,s )/(f i,max − f i,min ), i = {1, 2} and s = {1, 2,. .., P}, where f i,min and f i,max are the minimum and maximum value of the i th objective function; f i,s is the value of the objective function corresponding to solution s and objective i; P is the number of non-dominated solutions from the Pareto front.In order to demonstrate the merit of all the objectives corresponding to solution s, the normalized membership function µ * S is calculated [5]: The BCS corresponds to the maximum value of the index µ * S : µ max = Max(µ * S , s = 1, 2,. . ., P).

Classic SGO
The SGO is a metaheuristic, population-based algorithm that is inspired by human social group behavior and is used for solving complex problems [29].In the SGO, the population consists of a social group of N individuals X i = [x 1,i , x 2,i ,. .., x j,i ,. .., x n,i ] |i = 1,2,. .., N interacting and exchanging knowledge with each other, each individual representing a solution X i to the problem.The characteristics of the individuals represent the components x j,i , j = 1, 2,. .., n (n is the number of characteristics of an individual, which equals the size of the problem to be solved) of the solutions X i , and the ability of an individual to find a solution to the problem is measured by the fitness function f i = f (X i ).The SGO algorithm has three phases: initialization phase, improving phase, and acquiring phase.In the initialization phase, each component x j,i , j = 1, 2,. .., n of the solution X i , is randomly generated between the minimum (x min,j ) and maximum (x max,j ) limits: x j,i = x min,j + r 1 (x max,j − x min,j ) where r 1 is a random number uniformly distributed in the range (0, 1).
In the improving phase, each individual X i, i = 1, 2,. .., N of the group seeks to improve their traits x j,I by interacting with the best individual of the group at that moment, X best = [x 1 best , x 2 best ,. .., x j best ,. . .x n best ].Thus, the new traits (x j,i new ) of the individual X i new are determined with relation [29]: where c is the self-introspection parameter with values between 0 and 1, and r 2 is a uniformly generated random number in the range (0, 1).If the new vector X i new is better, then it is retained; otherwise, the old solution is kept.
In the acquiring phase, each individual X i aims to improve its traits (level of knowledge) through mutual, cumulative interactions with both the best individual X best and an individual randomly r (r ̸ = i) chosen from the population, X r = [x 1,r , x 2,r ,. .., x j,r ,. .., x n,r ].Depending on the quality of the individuals X i and X r , the new traits (x j,i new ) of the individual X i new are determined as follows [29] where, r 3 and r 4 denote independent arbitrary value in the (0, 1).
The solutions obtained using the SGO algorithm are iteratively improved using the updating relations from the improving and acquiring phases, until the maximum number of iterations t max is reached.A pseudo-code of the SGO algorithm is presented in detail in [29].

Modified SGO (MSGO)
The proposed MSGO algorithm includes the same phases as the original SGO algorithm, but some features of SGO have been modified to increase the efficiency of the MSGO algorithm.Two of these changes, presented below, aim at inserting two operators (one chaotic type and the other being highly disruptive polynomial (HDP)) in the structure of the solutions update relations to improve MSGO's ability to escape from local minima and obtain better quality solutions.
Chaotic operator: Chaos, due to its properties (unpredictability, non-periodic, ergodicity, non-converging, pseudo-randomness etc. [18]), has been successfully incorporated into the structure of various optimization algorithms (such as adaptive sparrow search algorithm [36] and moth-flame optimizer [37]) to cover a number of shortcomings related to low population diversity, slow convergence of the optimization process, stagnation in a local area etc.In this paper, chaos is included in the MSGO algorithm using a chaotic sequence (cx p ) generated by Logistic map [37] which can be effective in solving some EcD problems [18].The Logistic map is defined by relation [18]: where {cx p } p = 1. ..∞ represents the chaotic sequence generated by the Logistic map, at the p th iteration.The initial conditions are cx 0 ∈ (0, 1) and cx 0 ̸ = {0.0,0.25, 0.5, 0.75, 1}.HDP operator: In case of more complex mathematical functions, SGO cannot provide an efficient search in the whole solutions space by relying only on randomly generated sequences through numbers (r1-r4) or on randomly extracting a solution from the population (X r ).As a result, in order to increase the exploration and exploitation capacity of the MSGO algorithm, as well as the chances of exceeding local minima, a perturbation based on the HDP operator [38] is introduced into MSGO.Mathematically, the HDP operator is defined by relation [38]: where, ru is a random number uniformly generated in the range [0, 1]; LB i and UB i indicate the lower and upper boundary of the variables x j,i ; η is the index for the polynomial operator; α is a control parameter that is considered to have a value of 0.5 [39] (in this paper, the values of the parameter α are determined by experimental trials in the range (0, 1)); δ H is the equation that alternates moderate values with high values for δ, while δ L is the equation that generates low values for δ; δ 1 and δ 2 are calculated, based on the LB i and UB i imposed boundaries for the variable x j,i , with the relations In general, the HDP operator is used as a mutation operator that intervenes in altering some x j,i values of a X i solution with a certain predetermined probability.In MSGO, the HDP operator is not used as a mutation operator, but as an operator applied to modify each variable x j,i of a solution X i (practically, the mutation probability is considered equal to 1).Thus, depending on the HDP parameters (η, ru, δ 1 , α), the δ values generated by (26) may be high, which could cause some of the x j,i components of some X i solutions to reach the limits of the search space relatively quickly (the δ values could generate excessively large search steps).To avoid this problem, the δ values generated by ( 26) will be attenuated/smoothed according to the current iteration t, if δ exceeds a certain limit value (δ max ): Figure 1 shows the first 10,000 values of δ generated by relation (26), considering the following settings for the HDP operator parameters: The δ H values generated by ( 26) under the condition ru ≤ α = 0.5 are marked in blue, and the δ L values in red.Since the δ H values are much higher than δ L , in order to have a clearer picture, the δ L values have been detailed in a separate graph on the interval (1, 10,000).It is seen that the δ H equation alternates moderate values and high ones, while the δ L equation generates only low values.Figure 2 shows the attenuated δ a values obtained using (27).Initially, the δ a values are higher (being able to ensure an efficient search in the whole solutions space), then with the increase in the number of simulations they attenuate, the search step is reduced, making the transition from exploration to exploitation, and at the end of the process the values are equal to or lower than the δ max limit (which facilitates the exploitation phase).
Sustainability 2024, 16, x FOR PEER REVIEW 10 of 35 to 1).Thus, depending on the HDP parameters (η, ru, δ1, α), the δ values generated by ( 26) may be high, which could cause some of the xj,i components of some Xi solutions to reach the limits of the search space relatively quickly (the δ values could generate excessively large search steps).To avoid this problem, the δ values generated by ( 26) will be attenuated/smoothed according to the current iteration t, if δ exceeds a certain limit value (δmax): Figure 1 shows the first 10,000 values of δ generated by relation (26), considering the following settings for the HDP operator parameters: δ1 = rnd(1), δ2 = 1 − δ1; η = 4, α = 0.5, δmax = 1.5.The δH values generated by (26) under the condition ru ≤ α = 0.5 are marked in blue, and the δL values in red.Since the δH values are much higher than δL, in order to have a clearer picture, the δL values have been detailed in a separate graph on the interval (1, 10,000).It is seen that the δH equation alternates moderate values and high ones, while the δL equation generates only low values.Figure 2 shows the attenuated δa values obtained using (27).Initially, the δa values are higher (being able to ensure an efficient search in the whole solutions space), then with the increase in the number of simulations they attenuate, the search step is reduced, making the transition from exploration to exploitation, and at the end of the process the values are equal to or lower than the δmax limit (which facilitates the exploitation phase).The changes implemented in MSGO compared to SGO are presented below: to 1).Thus, depending on the HDP parameters (η, ru, δ1, α), the δ values generated by (26) may be high, which could cause some of the xj,i components of some Xi solutions to reach the limits of the search space relatively quickly (the δ values could generate excessively large search steps).To avoid this problem, the δ values generated by ( 26) will be attenuated/smoothed according to the current iteration t, if δ exceeds a certain limit value (δmax): Figure 1 shows the first 10,000 values of δ generated by relation (26), considering the following settings for the HDP operator parameters: δ1 = rnd(1), δ2 = 1 − δ1; η = 4, α = 0.5, δmax = 1.5.The δH values generated by (26) under the condition ru ≤ α = 0.5 are marked in blue, and the δL values in red.Since the δH values are much higher than δL, in order to have a clearer picture, the δL values have been detailed in a separate graph on the interval (1, 10,000).It is seen that the δH equation alternates moderate values and high ones, while the δL equation generates only low values.Figure 2 shows the attenuated δa values obtained using (27).Initially, the δa values are higher (being able to ensure an efficient search in the whole solutions space), then with the increase in the number of simulations they attenuate, the search step is reduced, making the transition from exploration to exploitation, and at the end of the process the values are equal to or lower than the δmax limit (which facilitates the exploitation phase).The changes implemented in MSGO compared to SGO are presented below: The changes implemented in MSGO compared to SGO are presented below: a.The improving phase of MSGO is similar to that of SGO but the sequences of numbers r 2 randomly generated in SGO by (23) are replaced by sequences of numbers generated by the attenuated DHP operator δ a by relation (27) in the form: b.The acquiring phase of MSGO is similar to that of SGO, except for the following three modifications: b1.The sequences of numbers r 3 randomly generated in SGO by (24a) are replaced by sequences of numbers generated by the attenuated DHP operator δ a by relation (27).b2.The sequences of numbers r 3 randomly generated in SGO by (24b) are replaced by chaotic sequences (cx) generated by the Logistic map with relation (25).b3.In SGO, switching between relations (24a) and (24b) is performed by considering the fitness of the competitive solutions X i and X r , based on the condition f (X i ) < f (X r ).In MSGO, this condition is replaced by a random one, having the form rnd(1) < β, where rnd( 1) is a uniformly generated random number in the range (0, 1); β is a value determined by experimental trials.
The three changes (b1)-(b3) made in MSGO are embodied in relations (29a) and (29b): If rnd(1) < β Then In Algorithm 1, the calculation steps for applying the MSGO algorithm are presented in detail, with the modifications (b1)-(b3) mentioned above.The main steps of the MSGO algorithm are summarized in a flowchart shown in Figure 3.
Find the best X best solution from the population; End For i Until t ≥ tmax {stopping criterion} The best solution X best and the fitness f(X best ) are retained.
The main steps of the MSGO algorithm are summarized in a flowchart shown in Figure 3.

Implementing the MSGO for the EcD, EmD, or EED Problems
This section shows how to implement the MSGO algorithm for solving EcD, EmD, or EED problems.The application of MSGO for solving EcD, EmD, or EED problems is similar, the difference being given by the type of objective function to be minimized: cost C(PT, PW) defined by (1) for EcD problems, emission E(PT, PW) defined by (16) for EmD problems or Ф(PT, PW) defined by (20) for EED problems.For the implementation of the MSGO algorithm, a solution i is considered to be represented by the vector PTWi = [PTW1,i, PTW2,i,..., PTWj,i,…, PTWNt+Nw,i]|i = 1, 2,… N which includes the powers of thermal units (PTWj,i|j = 1,2,... Nt) and wind units (PTWj,i|j = Nt+1, Nt+2,..., Nt+Nw).If the system contains only thermal units, then the vector PTWi consists only of the powers generated by these units.We also denote PTW best (t) as the best solution obtained using MSGO up to iteration t.

Stages of MSGO Implementation for the EcD, EmD, or EED Problem
The MSGO algorithm applied to solve the EcD, Em, or EED problems formulated in Section 2 includes the following steps:

Implementing the MSGO for the EcD, EmD, or EED Problems
This section shows how to implement the MSGO algorithm for solving EcD, EmD, or EED problems.The application of MSGO for solving EcD, EmD, or EED problems is similar, the difference being given by the type of objective function to be minimized: cost C(PT, PW) defined by (1) for EcD problems, emission E(PT, PW) defined by ( 16) for EmD problems or ϕ(PT, PW) defined by (20) for EED problems.For the implementation of the MSGO algorithm, a solution i is considered to be represented by the vector PTW i = [PTW 1,i, PTW 2,i, . .., PTW j,i, . .., PTW Nt+Nw,i ]| i = 1, 2,. . .N which includes the powers of thermal units (PTW j,i | j = 1,2,. . .Nt ) and wind units (PTW j,i | j = Nt+1, Nt+2,. .., Nt+Nw ).If the system contains only thermal units, then the vector PTW i consists only of the powers generated by these units.We also denote PTW best(t) as the best solution obtained using MSGO up to iteration t.

Stages of MSGO Implementation for the EcD, EmD, or EED Problem
The MSGO algorithm applied to solve the EcD, Em, or EED problems formulated in Section 2 includes the following steps: Step 1: Specify test system input data for the EcD, EmD, or EED problem: the number of thermal (Nt) and wind (Nw) units; cost coefficients for the thermal units (a, b, c, e, f ); parameters (c, k) to the Weibull distribution; the characteristic values of the wind speed (v in , v r , v out ); the rated output power of the wind units (PWr); active power limits for thermal and wind units (PT min , PT max ; PW min , PW max ); the cost coefficients (c d , c o , c u ) and emissions coefficients (e o , e u ) to the wind units; the B-loss coefficients (B ij , B 0i , B 00 ); load demand (P D ) and accuracy ε.
Step 2: Set the parameters of the MSGO algorithm: N and t max ; Step 3: Random initialization of a population of N solutions, represented by the vectors PTW i|i = 1,2,. . .N 3.1: t = 0; 3.2: Chaotic sequence cx are initialized using the logistic map; 3.3: Randomly generate an initial population with N solutions (PTW i (0) | i = 1,2,. .., N ) using relation (22).Each solution respects the constraints defined by relations ( 12)-( 14); The constraint ( 14) is handled by a heuristic procedure CHM presented in Section 4.2; 3.4: Evaluate the initial solutions PTW i (0) | i = 1,2,. .., N using relation ( 1), ( 16) or (20) associated with the EcD, EmD, or EED problems; 3.5: Find the best initial solution PTW best(0) and the objective function associated with the addressed problem (EcD, EmD, or EED).12) and ( 13): if the PTW j,i new(t) power is outside the limits, then CHM from Section 4.2 is applied; End For j; 4.7: Checking the equality constraint ( 14): the new solution PTW i new(t) is adjusted using the equality constraint handling mechanism from Section 4.2; 4.8: Evaluate the new solution PTW i new(t) using relation ( 1), ( 16) or (20) depending on the addressed problem (EcD, EmD or EED): If the new solution PTW i new(t) is better, then it is retained; otherwise, the old solution is maintained; 4.9: Update the best solution PTW best(t) ; End For i; Step 5: Update solutions PTW i (t) | i = 1,2,. .., N in the acquiring phase.

5.1:
For i = 1 To N Do 5.2: Randomly select a solution PTW r (t) , rϵ{1, 2,. .., N }, r ̸ = i; 5.3: If rnd(1) < β Then 5.4: For j = 1 to n Do 5.5: Generate a δ value of the HDP operator using (26); 5.6: Determine δ a using relation (27); 5.7: Updating PTW j,i (t) using (29a), getting the new components PTW j,i new(t) of the new solutions PTW i new(t) ; 5.8: Apply the procedure for handling inequality constraints ( 12) and ( 13) from Section 4.2; End For j; End If 5.9: Else 5.10: For j = 1 to n Do 5.11: Generate a chaotic value cx by (25); 5.12: Updating PTW j,i (t) by (29b), obtaining PTW j,i new(t) ; 5.13: Apply CHM from Section 4.2 for constraints ( 12) and ( 13); End For j; End Else 5.14: Checking the equality constraint ( 14): the new solution PTW i new(t) is adjusted using the equality constraint handling mechanism from Section 4.2; 5.15: Evaluate the new solution PTW i new(t) using relation ( 1), ( 16) or (20) depending on the addressed problem (EcD, EmD, or EED): if the new solution PTW i new(t) is better, then it is retained; otherwise the old solution is maintained; 5.16: Update the best solution PTW new(t) ; End For i; Step 6: Stop the process: the calculation process is stopped when the maximum number of iterations (t max ) is reached.{End For t} Step 7: Memorize the best solution: The best solution PTW best and the objective function associated with the problem addressed (EcD, EmD, or EED).

Constraints Handling Mechanism (CHM)
The CHM refers to the adjusting method of the solutions that do not satisfy the inequality and equality constraints of the EcD, EmD, or EED problems defined by relations ( 12)-( 14).In the case of inequality constraints, defined by relations ( 12) and ( 13), if the power of a thermal or wind unit is outside the imposed operating limits (PT min or PT max , and PW min or PW max , respectively), then the power of this unit is set with the exceeded limit [40].To handle the equality constraint, defined by ( 14), we apply a CHM procedure presented in [41].The CHM(PTW) procedure has a single parameter, represented by the power vector PTW, and finally it returns a feasible solution that respects relation (14) with the imposed error ε.The CHM(PTW) is called whenever a PTW vector is updated in the initializing, the improving, or the acquiring phases of the MSGO or SGO algorithms.

Case Studies
The efficiency of the MSGO algorithm was tested on medium-sized (10-unit) and large (40-unit) systems, having different characteristics, such as operating limits of the units, transmission line losses, valve-point effect, wind power.According to the characteristics of the systems, four cases were studied, described below.
Case 1 (C 1 ): This case analyses a 10-unit system taking into account the transmission line losses.The power demand is 2000 MW.The cost coefficients (a i , b i , c i , e i , f i , i = 1, 2,. .., 10), B-loss coefficients (B ij ), and emission coefficients (α i , β i , γ i , η i , δ i ) are taken from [42].Case 2 (C 2 ): The second case is a system having 40 units that considers the valve-point effects.The system is analyzed without considering transmission line losses.The power demand is 10,500 MW.The cost coefficients are considered from [33] and emission coefficients are provided by [43].Case 3 (C 3 ): The case C 3 is a system with 40 units similar to case C 2 (the cost and emission coefficients are identical to those in case C 2 [33,43]), but transmission losses are considered, as well.The power demand is 10,500 MW.The B-loss coefficients are taken from [44].Case 4 (C 4 ): A 40-unit system derived from case C 3 by replacing the first two thermal units (PT 1 and PT 2 ) with the two wind units (PW 1 and PW 2 ).Each wind unit has a nominal power of 550 MW, and the minimum and maximum capacities are PW min,1 = PW min,2 = 0, PW max,1 = PW max,2 = 550 MW.The power demand is P D = 10,500 MW.We consider that the wind speed has a Weibull distribution.The shape and scale parameters (k, c) corresponding to the sites of the two wind units have the values [13]: k 1 = 1.5; c 1 = 15; k 2 = 1.5; c 2 = 15.Other wind-related characteristics have the following values [13]: cut-in wind speed (v in = 5 m/s), rated wind speed (v r = 15 m/s), cut-out wind speed (v out = 45 m/s), the cost coefficients associated with underestimation (c u 1 = c u 2 = 5) and overestimation (c o 1 = c o 2 = 5).The cost and emission coefficients for the thermal units (P 3 -P 40 ) and the B-loss coefficients are identical to those mentioned in C 3 case, being taken from [44].The characteristics of the thermal units and the B-loss coefficients for the 10-unit and 40-unit systems are presented in Appendix A, Tables A2-A5.
Each of the cases mentioned in Table 1 (C 1 -C 4 ) were studied based on the mathematical optimization models (presented in Section 2) of the problems EcD (cases C 1a -C 4a ), EmD (cases C 1b -C 4b ), and EED (cases C 1c -C 4c ).The SGO and MSGO algorithms perform 50 independent runs for each case study, retaining four statistical indicators: Best cost/emission (B), Average cost/emission (A), Worst cost/emission (W), and the standard deviation (SD).Algorithms SGO and MSGO have been implemented in MathCAD and run on a PC with an Intel i5 processor, 2.2 GHz CPU and 4 GB of RAM.

Setting the Parameters
Adjusting the parameters of metaheuristic algorithms is an essential task with multiple positive effects, such as improved performance, adaptation of the algorithms to a specific problem, better execution times, and stability, etc.The proposed MSGO algorithm has three specific parameters (α, β, and δ max ).A procedure based on the design of an experimental plan is used to set the specific parameters.They are successively set (for example, in the following order α, β, and δ max ), starting from a set of initial values [45].Thus, one parameter is varied, and the others are fixed either with the initial values or with the already set values.This procedure allows the evaluation of the MSGO algorithm performance through testing different combinations of the interest parameter values, followed by the selection of the best combination based on a predetermined criterion.The test values considered for setting the specific parameters are α = {0, 0.25, 0.5, 0.75, 1}, β = {0, 0.25, 0.5, 0.75, 1}, and δ max = {1.2,1.5, 2}.The initial values are α = 0, β = 0, and δ max = 1.2.The parameter selection criterion is based on Average cost/emission obtained from 25 runs for each case study.Note that the experimental analysis does not guarantee the best values for the parameters of MSGO algorithm.However, the results obtained using MSGO following this selection process show that the algorithm parameters were reasonably set Table 1 shows the parameters of the MSGO algorithm (N, t max , α, β, δ max ) for all the analyzed cases.In the case of the SGO algorithm, the values of the common parameters (N and t max ) are the same as the values set for the MSGO algorithm in each case analyzed (these are shown in Table 1).The SGO algorithm has only one specific parameter c, the value of which being set to the one recommended in [29] (c = 0.2).For the fair comparison of the SGO and MSGO algorithms, the number of evaluations of the objective functions (NE) is considered equal in all analyzed cases (C 1 -C 4 ), which is mentioned in Table 1.

Results for EcD (Cases C 1a -C 4a ) and EmD (Cases C 1b -C 4b ) Problems
The best solutions obtained using SGO and MSGO: The optimal scheduling of the thermal units for the cases C 1a , C 1b (10 units), and C 2a , C 2b (40 units) obtained using the MSGO and SGO algorithms are presented in Tables 2 and 3, respectively.
In the cases C 1a and C 1b , the system being of relatively small sizes, the MSGO and SGO algorithms find approximately the same optimal operating solution, and as a result, the derived quantities (Cost, Emission, P L ) will be approximately the same.However, mathematically, MSGO performs slightly better than SGO.Thus, in the case of C 1a the best costs obtained using MSGO and SGO are Cost C1a(MSGO) = 111,497.6301$/h and Cost C1a(SGO) = 111,497.6302$/h, respectively.For the C 1b case, the best emissions found by MSGO and SGO are Emission C1b(MSGO) = 3932.2432lb/h and Emission C1b(SGO) = 3932.2433lb/h, respectively.Also, in the case of C 2b (    On the contrary, in the case of C 2a (where the functions that model the cost have a higher complexity and determine a larger number of local minima), the MSGO algorithm can identify a much better solution than SGO (Cost C2a(MSGO) = 121,426.7039$/h < Cost C2a(SGO) = 121,509.8209$/h).
Table 4 shows the best solutions in terms of cost and emissions obtained using the MSGO algorithm for the cases C 3a , C 3b (without wind), C 4a , and C 4b (with wind).The solutions for C 4a and C 4b cases indicate that the two wind units (PW 1 , PW 2 ) are scheduled to operate at full capacity (550 MW) in both the EcD and EmD problems.Thus, according to the statement in case C 4 , the wind units will replace two thermal units PT 1 , PT 2 (that operate at maximum capacity PT 1 ≈ PT 2 ≈ 114 MW in C 3a and C 3b cases), and will also reduce the operating power level of other thermal units (PT 3 -PT 40 ).Based on the EcD and EmD models, it can be seen that if the wind units operate at full capacity it results in an E u power very close to zero (E u 1 (PW 1 = 550) ≈ E u 2 (PW 2 ) ≈ 0 MW), so their corresponding costs and emissions will be close to zero.In contrast, the average powers associated with overestimation will tend towards maximum values 7635 MW , while their corresponding costs and emissions are 2 × 1153.817$/h, and 2 × 1153.817t/h.The risks taken by the system operator when planning the wind units to operate at values close to the maximum capacity are covered by the reduction of costs and emissions from the thermal units (due to the reduction of the PT 3 -PT 40 powers in C 4a (C 4b ) cases compared to C 3a (C 3b )).
Comparing the solutions for the cases with wind and those without wind, it can be seen that the inclusion of the wind units (PW 1 and PW 2 ) results in a substantial reduction in cost (from Cost C3a(MSGO) = 136,454.33693$/h to Cost C4a(MSGO) = 123,161.88666$/h, representing 9.74%), and emissions (from Emission C3b(MSGO) = 347,578.49052t/h to Emission C4b(MSGO) = 193,311.54070t/h representing 44.38%).Also, power losses are reduced from 958.62 MW in the case of C 3a to 936.70 MW in the case of C 4a (2.28%) and from 1040.54 MW in the case of C 3b to 1009.74 MW in the case of C 4b (2.96%), respectively.
Comparison of MSGO algorithm with other algorithms: In Tables 5-12 the values of the statistical items (B, A, W, and SD) are presented.They are obtained using different algorithms for EcD problem (cases C 1a -C 4a ) and EmD problem (C 1b -C 4b ), respectively.
Analyzing the values from the mentioned tables, it can be seen that in all studied cases the MSGO algorithm obtains statistical indicators as good or better than the competing algorithms mentioned in these tables.The exception is the Worst cost item for the algorithms Jaya-SML [46] (case C 2a -Table 6), ORCCRO [47] (case C 3a -Table 7), DE (case C 4a -Table 8), and the SD item for the algorithms CSS [48] (cases C 3a -Table 7), DE, and SCA (cases C 4a -Table 8).
We specify that the specific parameters of the PSO, DE, and SCA algorithms were set by performing several experiments in which the number of evaluations is considered to be the same as for the SGO and MSGO algorithms (shown in Table 1).Thus, the settings performed for these algorithms applied in the cases C 3b , C 4a , and C 4b are the following:

•
For the PSO algorithm (N = 50, t max = 400, c 1 = c 2 = 2, w min = 0.3, w max = 0.9 in cases C 3b and C 4b , respectively, N = 50, t max = 700, c 1 = c 2 = 2, w min = 0.3, w max = 0.9 in the case of C 4a ; where c 1 and c 2 are acceleration coefficients, w min and w max are the initial and final inertial weights).

•
For the DE algorithm (N = 50, t max = 400, CR = 0.2, F = 0.4 in cases C 3b and C 4b , respectively, N = 50, t max = 700, CR = 0.1, F = 0.8 in the case of C 4a ; where, CR is the crossover rate, and F is the scaling factor).
The comparison of the MSGO and SGO algorithms: Tables 5-12 present the values of the statistical items (B, A, W, SD) obtained using the SGO and MSGO algorithms for all the study cases (C 1a -C 4a and C 1b -C 4b , respectively).In all situations, the statistical items obtained using MSGO are higher than those obtained using SGO (except for the cases mentioned in the previous paragraph), but there are also cases where the differences between the SGO and MSGO algorithms are small (in particular, the cases C 1b -C 4b associated with EmD problem, to which case C 1a can be added).For this reason, the SGO and MSGO algorithms are compared using the non-parametric Wilcoxon statistical test, considering a significance-level of 1%.Also, 50 values/algorithm were simulated to compare SGO and MSGO, resulting in 50 pairs of values that are compared using the Wilcoxon test.
Table 13 shows the statistical item p-value after applying the Wilcoxon test for the comparison of SGO and MSGO in all cases studied.For the EcD (C 1a -C 4a ) and EmD (C 1b -C 4b ) problems, the p-value is less than 0.01 (except for cases C 2b and C 3b ), which indicates that MSGO is statistically significantly better than the SGO algorithm (getting six wins out of a possible eight).In order to investigate the ability of MSGO to identify a better solution compared to SGO, we will evaluate the cost savings, as well as the emission reduction for cases C 1a -C 4a and C 1b -C 4b .Thus, for C 1b -C 4b (Tables 9-12) the MSGO and SGO algorithms perform very well at solving the EmD problem.As a result, the emission reduction achieved by MSGO compared to SGO is insignificant for this type of problem.The situation is similar in the EcD problem, case C 1a (Table 5), when both algorithms perform very well on the small 10-unit system.However, the cost savings and emission reduction may be important when comparing MSGO with other competing algorithms.Thus, for C 1a case (Table 5), analyzing the Average cost item, the average cost reduction is 10,124.12$/h compared to QPSO [59].Also, in the EmD problem (Tables 9-12), based on the Average emission item, the average hourly emission reduction is 109.6 lb/h (compared to QPSO [59], in case C 1b ), 542.3 t/h (relative to DE, in case C 3b ), and 1220.5 t/h (relative to DE, in C 4b ).Following the results in Tables 6-8, the cost savings achieved by MSGO compared to SGO is significant: 368.16 $/h (for case C 2a ), 416.80 $/h (case C 3a ), and 525.04 $/h (case C 4a ).These cost reductions are maintained or increased when comparing MSGO with other algorithms presented in Tables 6-8.
Convergence process: Figures 4 and 5 show the convergence characteristics in terms of cost and emissions obtained using the MSGO algorithm for all studied cases C 1a -C 4a and C 1b -C 4b .In all analyzed situations, convergence is fast in the first iterations (approximately 15% of the maximum number of iterations), and in the last iterations the convergence is slow, the iterative process being stabilized.
comparing MSGO with other competing algorithms.Thus, for C1a case (Table 5), analyzing the Average cost item, the average cost reduction is 10,124.12$/h compared to QPSO [59].Also, in the EmD problem (Tables 9-12), based on the Average emission item, the average hourly emission reduction is 109.6 lb/h (compared to QPSO [59], in case C1b), 542.3 t/h (relative to DE, in case C3b), and 1220.5 t/h (relative to DE, in C4b).Following the results in Tables 6-8, the cost savings achieved by MSGO compared to SGO is significant: 368.16 $/h (for case C2a), 416.80 $/h (case C3a), and 525.04 $/h (case C4a).These cost reductions are maintained or increased when comparing MSGO with other algorithms presented in Tables 6-8.
Convergence process: Figures 4 and 5 show the convergence characteristics in terms of cost and emissions obtained using the MSGO algorithm for all studied cases C1a-C4a and C1b-C4b.In all analyzed situations, convergence is fast in the first iterations (approximately 15% of the maximum number of iterations), and in the last iterations the convergence is slow, the iterative process being stabilized.

Chart Title
Fuel cost (C1a)   In the case of emission convergence characteristics (C1b-C4b), the iterative process stabilizes much faster than in the case of cost characteristics (C1a-C4a).Thus, in the case of the EmD problem (C1b-C4b), the quasi-stable process starts after about 20-30% of the maximum number of iterations, and in the EcD problem after approximately 40-50% of the maximum number of iterations.We mention that the stable values towards which the convergence characteristics of the MSGO algorithm tend are mentioned for each case in Tables 5-12.
The accuracy and computation time: The accuracy (∆P) of the computations refers to the inequality and equality constraints of the optimization model presented in Section 2 and can be seen for the best solutions presented in Tables 2-4.It is observed that the constraints are respected, the size of ∆P being less than 10 -5 MW for all cases studied.Also, for each case (C1a-C4a and C1b-C4b), the average execution time (Time) obtained using the MSGO algorithm is indicated in Tables 2-4.In cases C1a and C1b, the execution time is very good (under 0.2 s); for cases C2a and C2b it is higher (under 10 s) due to the larger size of the analyzed system (40-unit).In the cases C3a-C4a and C3b-C4b, due to the large size of the analyzed systems (40-unit) and the need to calculate power losses, the average execution time goes up to 33 s.In the case of emission convergence characteristics (C 1b -C 4b ), the iterative process stabilizes much faster than in the case of cost characteristics (C 1a -C 4a ).Thus, in the case of the EmD problem (C 1b -C 4b ), the quasi-stable process starts after about 20-30% of the maximum number of iterations, and in the EcD problem after approximately 40-50% of the maximum number of iterations.We mention that the stable values towards which the convergence characteristics of the MSGO algorithm tend are mentioned for each case in Tables 5-12.
The accuracy and computation time: The accuracy (∆P) of the computations refers to the inequality and equality constraints of the optimization model presented in Section 2 and can be seen for the best solutions presented in Tables 2-4.It is observed that the constraints are respected, the size of ∆P being less than 10 −5 MW for all cases studied.Also, for each case (C 1a -C 4a and C 1b -C 4b ), the average execution time (Time) obtained using the MSGO algorithm is indicated in Tables 2-4.In cases C 1a and C 1b , the execution time is very good (under 0.2 s); for cases C 2a and C 2b it is higher (under 10 s) due to the larger size of the analyzed system (40-unit).In the cases C 3a -C 4a and C 3b -C 4b , due to the large size of the analyzed systems (40-unit) and the need to calculate power losses, the average execution time goes up to 33 s.

Results for EED Problem (Cases C 1c -C 4c )
The EED problem considers the cases C 1c -C 4c , in which the single-objective function ϕ defined by relation (20) is minimized.In these cases, it is of interest to determine the Pareto front solutions, as well as the best compromise solution (BCS) obtained using the SGO and MSGO algorithms.
To estimate the Pareto front, the weighting factor ω in relation ( 20) is varied between 0 and 1 with a step of 0.1 (in the end 11 points are obtained in the Cost-Emission objective plan).However, in order to obtain a Pareto front as uniform as possible and to show graphically that the solutions obtained using MSGO dominate the solutions obtained using other algorithms, another 12 points have been added to the 11 points for the case C 1c (the added points correspond to the values ω = {0.25,0.35, 0.41, 0.42, 0.43, 0.435, 0.436, 0.44, 0.441, 0.445, 0.47, 0.55}) and 3 another points in the case of C 2c (ω = {0.51,0.52, 0.55}), respectively.In the case of C 3c and C 4c , the number of points considered is 11.
Figures 6-8 show the Pareto fronts obtained using the SGO and MSGO algorithms for the cases C 1c -C 4c .Also, the best compromise solutions obtained using SGO, MSGO, and other competing algorithms applied in the literature are indicated.
the best compromise solution determined by SGO and MSGO.Comparison of solutions from the Pareto front of MSGO with BCSs obtained using other algorithms: In the case of C1c, in order to test the ability of MSGO to identify solutions from the Pareto front, it was compared with BCSs obtained using a group of algorithms (denoted G1) presented in the literature G1 = {EMOCA [24], MODE [42], NSGA II [42], TLBO [26], GSA [70], PDE [42], SMA [71], εv-MOGA [49], SPES-2 [42]}.To highlight the Pareto front solutions obtained using MSGO and BCSs obtained using competing algorithms in the G1 group, in Figure 6, a detail of the Cost-Emission plan is shown.In this detail of Figure 6, it can be visually observed that the Pareto front solutions obtained using MSGO (marked with red circles) dominate, in terms of Cost and Emission, all BCSs obtained using competing algorithms G1 (marked with blue circles).Thus, the points marked with red circles (obtained using MSGO) will remove from the Pareto front of MSGO all the points marked by blue circles (obtained using competing G1 algorithms).The exception is the IKSO algorithm [72], which obtains a non-dominated solution using the solutions in the Pareto front of the MSGO (a discussion on IKSO [72] will be made below).In case of C2c, in Figure 7, a detail of the Cost-Emission plan is shown highlighting the Pareto front solutions obtained using MSGO, and BCSs obtained through several algorithms presented in the literature, denoted G2, G2 = {LFA [16], MODE [42], SPEA-2 [42], εv-MOGA [49], PDE [42], NSGA II [22], TLBO [26], CSOA [20], KKO [21], KSO [22]  Best compromise solution by SGO and MSGO: The BCSs obtained using the SGO and MSGO algorithms are identified using a fuzzy-based mechanism presented in Section 2.3.The results obtained through this mechanism indicate that BCS corresponds to a weighting factor equal to ω = 0.5, both for SGO and MSGO.
As a result, in Table 14, we present for each case C 1c -C 4c , the values for Cost, Emission, and item µ max (determined according to the methodology in Section 2.3) corresponding to the best compromise solution determined by SGO and MSGO.
Comparison of solutions from the Pareto front of MSGO with BCSs obtained using other algorithms: In the case of C 1c , in order to test the ability of MSGO to identify solutions from the Pareto front, it was compared with BCSs obtained using a group of algorithms (denoted G1) presented in the literature G1 = {EMOCA [24], MODE [42], NSGA II [42], TLBO [26], GSA [70], PDE [42], SMA [71], εv-MOGA [49], SPES-2 [42]}.To highlight the Pareto front solutions obtained using MSGO and BCSs obtained using competing algorithms in the G1 group, in Figure 6, a detail of the Cost-Emission plan is shown.In this detail of Figure 6, it can be visually observed that the Pareto front solutions obtained using MSGO (marked with red circles) dominate, in terms of Cost and Emission, all BCSs obtained using competing algorithms G1 (marked with blue circles).Thus, the points marked with red circles (obtained using MSGO) will remove from the Pareto front of MSGO all the points marked by blue circles (obtained using competing G1 algorithms).The exception is the IKSO algorithm [72], which obtains a non-dominated solution using the solutions in the Pareto front of the MSGO (a discussion on IKSO [72] will be made below).the image it can be visually observed that the Pareto front solutions obtained using MSGO (marked with red circles) dominate, in terms of Cost and Emission, the BCSs obtained using the G2 group algorithms (marked with blue circles).Comparison of BCSs obtained using MSGO and SGO: To compare BCSs obtained using the SGO and MSGO algorithms, we will include in the Pareto front of MSGO the best compromise solution identified for the SGO algorithm.Similarly, in the case of C1c, the BCS identified using the IKSO algorithm [72] will be included in the Pareto front of MSGO.Note that BCSs obtained using the SGO algorithm (cases C1c-C4c) or IKSO [72] (in the case of C1c) must be non-dominated solutions in the MSGO Pareto front.Thus, MSGO will have a Pareto front increased by one value (in cases C3c, corresponding to the BCS of SGO), respectively by two values (in cases C1c, corresponding to the BCSs of SGO and IKSO [72]).In cases C2c and C4c, the BCS of MSGO dominates the BCS of SGO in terms of Cost and Emission (see Table 14).The values from these extended Pareto fronts (with one or two other values) are subject to the methodology for establishing the BCS presented in Section 2.3, and the results of interest are given in Table 15.Table 15 shows the values of the item µmax obtained using the SGO, MSGO, and IKSO algorithms (we specify that the IKSO algorithm [72] is used only in the case of C1c).From Table 15, it can be seen that, for the cases C1c and C3c, the maximum value of the item µ corresponds to the MSGO the image it can be visually observed that the Pareto front solutions obtained using MSGO (marked with red circles) dominate, in terms of Cost and Emission, the BCSs obtained using the G2 group algorithms (marked with blue circles).Comparison of BCSs obtained using MSGO and SGO: To compare BCSs obtained using the SGO and MSGO algorithms, we will include in the Pareto front of MSGO the best compromise solution identified for the SGO algorithm.Similarly, in the case of C1c, the BCS identified using the IKSO algorithm [72] will be included in the Pareto front of MSGO.Note that BCSs obtained using the SGO algorithm (cases C1c-C4c) or IKSO [72] (in the case of C1c) must be non-dominated solutions in the MSGO Pareto front.Thus, MSGO will have a Pareto front increased by one value (in cases C3c, corresponding to the BCS of SGO), respectively by two values (in cases C1c, corresponding to the BCSs of SGO and IKSO [72]).In cases C2c and C4c, the BCS of MSGO dominates the BCS of SGO in terms of Cost and Emission (see Table 14).The values from these extended Pareto fronts (with one or two other values) are subject to the methodology for establishing the BCS presented in Section 2.3, and the results of interest are given in Table 15.Table 15 shows the values of the item µmax obtained using the SGO, MSGO, and IKSO algorithms (we specify that the IKSO algorithm [72] is used only in the case of C1c).From Table 15, it can be seen that, for the cases C1c and C3c, the maximum value of the item µ corresponds to the MSGO   In case of C 2c , in Figure 7, a detail of the Cost-Emission plan is shown highlighting the Pareto front solutions obtained using MSGO, and BCSs obtained through several algorithms presented in the literature, denoted G2, G2 = {LFA [16], MODE [42], SPEA-2 [42], εv-MOGA [49], PDE [42], NSGA II [22], TLBO [26], CSOA [20], KKO [21], KSO [22]}.From the image it can be visually observed that the Pareto front solutions obtained using MSGO (marked with red circles) dominate, in terms of Cost and Emission, the BCSs obtained using the G2 group algorithms (marked with blue circles).
Comparison of BCSs obtained using MSGO and SGO: To compare BCSs obtained using the SGO and MSGO algorithms, we will include in the Pareto front of MSGO the best compromise solution identified for the SGO algorithm.Similarly, in the case of C 1c , the BCS identified using the IKSO algorithm [72] will be included in the Pareto front of MSGO.Note that BCSs obtained using the SGO algorithm (cases C 1c -C 4c ) or IKSO [72] (in the case of C 1c ) must be non-dominated solutions in the MSGO Pareto front.Thus, MSGO will have a Pareto front increased by one value (in cases C 3c , corresponding to the BCS of SGO), respectively by two values (in cases C 1c , corresponding to the BCSs of SGO and IKSO [72]).In cases C 2c and C 4c , the BCS of MSGO dominates the BCS of SGO in terms of Cost and Emission (see Table 14).The values from these extended Pareto fronts (with one or two other values) are subject to the methodology for establishing the BCS presented in Section 2.3, and the results of interest are given in Table 15.Table 15 shows the values of the item µ max obtained using the SGO, MSGO, and IKSO algorithms (we specify that the IKSO algorithm [72] is used only in the case of C 1c ).From Table 15, it can be seen that, for the cases C 1c and C 3c , the maximum value of the item µ corresponds to the MSGO algorithm (the weighting factor being for each case ω = 0.5).As a result, in all cases (C 1c -C 4c ), the MSGO algorithm can identify a better compromise solution than SGO.Figures 5-8 show the Pareto fronts obtained using SGO and MSGO and indicate the BCSs identified using these algorithms.-the IKSO algorithm [72] is not used in the C 2c -C 4c cases.The bold values correspond to the winning algorithm.
Comparison of SGO and MSGO Pareto fronts: Two metrics are used to evaluate and compare the quality of the Pareto fronts obtained using the SGO and MSGO algorithms: C-metric [73] and hyper-volume [74].The C-metric, denoted C(S 1 , S 2 ) and applied to the non-dominated solution sets S 1 and S 2 , indicates the percentage of solutions from the S 2 set dominated by the solutions from S 1 .For example, if C(S 1 , S 2 ) = 100%, all solutions in S 2 are dominated by those in S 1 .Currently, C(S 1 , S 2 ) ̸ = C(S 2 , S 1 ) is required to calculate both metrics C(S 1 , S 2 ) and C(S 2 , S 1 ).The hyper-volume metric measures the volume (area in the case of two objectives) that is dominated by the Pareto front and is located below a predetermined reference point.Typically, higher values of the hyper-volume metric are associated with a better-quality Pareto front.
Table 16 shows, for each case, the values of the C-metric and hyper-volume metrics corresponding to the solution sets obtained using the SGO and MSGO algorithms.From Table 16 it can be seen that the C-metric (SGO, MSGO) = 0, which indicates that the solutions in the Pareto front obtained using SGO do not dominate any solution in the Pareto front obtained using the MSGO algorithm.Instead, the inverse C-metric (MSGO, SGO) indicates the percentages of SGO solutions dominated by MSGO solutions.Also, the hyper-volume metric shows higher values in the case of MSGO compared to SGO.Considering the values obtained using the two metrics, we estimate that MSGO can obtain a Pareto front of better quality than SGO.Average execution time: Tables 2-4 show the values of the average execution time (Time) obtained using the MSGO algorithm for each case (C 1c -C 4c ) of the EED problem.It can be seen that the values obtained for the EED problem are slightly higher than for the EcD or EmD problems.

Conclusions
In this article, a modified SGO is proposed in which several terms in the solution update relations are disturbed by including chaotic sequences generated by the Logistic map and sequences generated by the HDP operator.Also, switching between solution update relations (in the acquiring phase) is achieved by introducing a random condition instead of a condition based on the value of the fitness function of the competing solutions.MSGO has been successfully tested for solving EcD, EmD, and EED problems for medium and large power systems.
The effectiveness of MSGO was tested on four cases for each type of problem (EcD, EmD, and EED).For all case studies targeting EcD, and EmD problems, MSGO obtained solutions of better or equal quality than well-known algorithms (such as: DE, ABC, PSO, SCA, etc.) or improved varieties (except for the cases mentioned in Section 5.2).Also, the convergence process was fast and the stability of the MSGO algorithm evaluated by the SD item was very good in cases C 1b -C 4b and relatively good in the cases C 1a -C 4a .The statistical items (B, A, W and SD) and the results of the Wilcoxon test show that MSGO is more efficient than SGO in six cases (out of a possible eight) and equally good in two cases (C 2b and C 3b in Table 13), indicating that the changes made in MSGO had a positive effect.In the case of the EED problem, MSGO was able to extract a better compromise solution than SGO or other algorithms for all studied cases C 1c -C 4c .Moreover, the values of the hyper-volume and C-metric indicate that MSGO can obtain a Pareto front of superior quality compared to SGO.It should be mentioned that the inclusion of wind sources (in the case of C 4a /C 4b compared to C 3a /C 3b ) brings benefits both in terms of cost reduction (about 10%) and polluting emissions (about 45%).The average computation time of MSGO and SGO algorithms is similar, and the efficiency of both algorithms is good in all case studies.By applying economic and/or emissions dispatching using MSGO, electricity producers can make informed decisions to operate a sustainable energy system.

Figure 1 .
Figure 1.HDP values generated for δ H and δ L .

Figure 2 .
Figure 2. The attenuated values obtained for δ a .

Figure 6 .
Figure 6.Pareto fronts and BCSs obtained using SGO, MSGO, and other algorithms for case C 1c .

Figure 7 .
Figure 7. Pareto fronts and BCSs obtained using SGO and MSGO for case C2c.

Figure 7 .
Figure 7. Pareto fronts and BCSs obtained using SGO and MSGO for case C2c.

Figure 8 .
Figure 8. Pareto fronts and BCSs obtained using SGO and MSGO: (a) for case C 3c ; (b) for case C 4c .

Author Contributions:
Conceptualization, D.C.S., C.H. and M.L.S.; software data curation, D.C.S.; investigation, D.C.S., C.H., M.L.S., G.B., C.B. and F.C.D.; writing-original draft preparation, D.C.S., M.L.S. and C.B.; writing-review and editing, C.H., G.B., C.B. and F.C.D.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by the University of Oradea within the framework of the grant competition "Scientific Research of Excellence Related to Priority Areas with Capitalization through Technology Transfer: INO-TRANSFER-UO-2nd Edition" projects no.239/2022 and no.236/2022.

Table 1 .
Characteristics of the analyzed systems and values set for the SGO and MSGO parameters.

Table 2 .
Generation schedule obtained using SGO and MSGO for cases C 1a , C 1b, and C 1c (P D = 2000 MW).
* ∆P represents the accuracy with which the power balance is satisfied: ∆P = ΣPT i − P L − P D .

Table 3 .
Generation schedule obtained using SGO and MSGO for cases C 2a , C 2b, and C 2c (P D = 10,500 MW).

Table 3 .
Cont. ∆P represents the accuracy with which the power balance is satisfied: ∆P = ΣPT i − P L − P D .

Table 4 .
Generation schedule obtained using MSGO for cases C 3a , C 3b, and C 3c (with loss of power and without wind) and C 4a , C 4b , C 4c (with loss of power and wind).

Case C 4a (With Wind) Case C 3b (Without Wind) Case C 4b (With Wind) Case C 3c (Without Wind) Case C 4c (With Wind)
∆P represents the accuracy with which the power balance is satisfied: ∆P = (ΣPT i + ΣPT i ) − P L − P D .

Table 5 .
Values for the statistical items obtained using different algorithms for the cases C 1a -Best Cost (10 units with losses, 2000 MW).
Indicates the cost saving through MSGO compared to other algorithms using the Average Cost item. *

Table 6 .
Values of the statistical items obtained using different algorithms for the cases C 2a -Best Cost (40 units without losses, 10,500 MW).

Table 7 .
Values of the statistical items obtained using different algorithms for the cases C 3a -Best Cost (40 units with losses and without wind, 10,500 MW).

Table 8 .
Values of the statistical items obtained using different algorithms for the cases C 4a -Best Cost (40 units with losses and with wind, 10,500 MW).

Table 9 .
Values of the statistical items obtained using different algorithms for the cases C 1b -Best Emission (10 units with losses, 2000 MW).
* Indicates the emisssion reduction through MSGO compared to other algorithms using the Average Cost item; ** in the case of the DE algorithm, the correct value is 3932.41728lb/h.

Table 10 .
Values of the statistical items obtained using different algorithms for the cases C 2b -Best Emission (40 units without losses, 10,500 MW).

Table 11 .
Values of the statistical items obtained using different algorithms for the cases C 3b -Best Issue (40 units with losses and without wind, 10,500 MW).

Table 12 .
Values of the statistical items obtained using different algorithms for the cases C 4b -Best Issue (40 units with losses and with wind, 10,500 MW).

Table 13 .
Comparison of SGO and MSGO by the Wilcoxon test.− and R + represent the mean of the negative and positive ranks, respectively. R

Table 14 .
The values corresponding to cost and emission obtained using SGO and MSGO for the cases C1c-C4c.BCS dominates SGO's BCS in terms of Cost and Emission; ** MSGO Cost savings and emission reduction comparing to SGO, for each case study BSCs.
}. From Pareto fronts and BCSs obtained using SGO and MSGO for case C 2c .

Table 14 .
The values corresponding to cost and emission obtained using SGO and MSGO for the cases C 1c -C 4c .
* MSGO's BCS dominates SGO's BCS in terms of Cost and Emission; ** MSGO Cost savings and emission reduction comparing to SGO, for each case study BSCs.

Table 15 .
Values of item µ max obtained using SGO, MSGO and IKSO for cases C 1c -C 4c .

Table 16 .
Values of the metrics for assessing the Pareto fronts of SGO and MSGO.When analyzingTable 4 from the EED problem point of view-case C 4c -it shows that the BCS is obtained when the wind units operate at maximum capacity (PW 1 ≈ PW 2 ≈ 550 MW).When comparing the BCS-including wind (case C 4c ) with the BCS-without wind (case C 3c ), it results in lower costs (by 8.9%) and emissions (by 38.3%) for case C 4c .Table 14 shows Cost savings and Emission reduction by MSGO compared to SGO for BCSs obtained in cases C 1c -C 4c .Positive values indicate that MSGO achieves better cost or lower emissions than SGO, and negative values indicate the opposite situation.

Table A2 .
Thermal units' characteristics for the 10-unit test system.

Table A4 .
Thermal units' characteristics for the 40-unit test system.