A High-Performance Stochastic Fractal Search Algorithm for Optimal Generation Dispatch Problem

This paper proposes applications of a modified stochastic fractal search algorithm (MSFS) to solve the economic load dispatch problem (ELD) in which valve-point effects, prohibited operating zones, power losses in all conductors, multi-fuel sources and other constraints of power system are taken into consideration. The proposed method is first developed in the study by performing two modifications on two procedures of new solution generation from conventional stochastic fractal search (SFS). The first modification is used to change the strategy of producing new solutions of the first and the second update procedures while the second one is to newly update the worst solutions in the first update process and the best solutions in the second update process. These modifications have major influence on the solution search performance of the proposed method. All improvements of the proposed method can be illustrated by solving and analyzing results from various test systems with different system scales including 3-unit, 6-unit, 10-unit, and 20-unit systems. Comparison of results obtained by MSFS, SFS, and other existing methods indicates that the proposed MSFS method is more effective and robust than compared methods in terms of solution quality, high-quality solution search stability and convergence process. Consequently, the proposed method should be used as a very favorable optimization method for the ELD problem and it should be tried for other optimization problems in electrical engineering.


Introduction
One of the most conspicuous trends in the 21st century is the economic growth, leading to ever-growing needs of energy consumption in all of the activities including production, services, and daily domestic households.They necessitate a power system to supply an adequate amount of energy in order to address rising power demands.However, one of the most complicated issues of power system supply is suitable electricity energy management that can use fuel cost-effectively and satisfy the system constraints exactly.The objective and constraints of such power system are taken into account in the economic load dispatch problem (ELD) [1].
Typically, the simplest cost-power feature of each thermal unit in the ELD problem is approximately represented as a quadratic function since valve-point effects and prohibited operation zones are not (SFS), which is also a population-based meta-heuristic algorithm and was derived from the natural evolution process [43].SFS has received much interest from researchers and it has been successfully utilized to resolve optimization problems in engineering fields so far.The main goal of this method was to create an efficient solution search strategy by using diffusion procedure and two update mechanisms.This algorithm can find near-optimal solutions and overtake the disadvantages from other meta-heuristic methods, such as premature convergence in a local optimum zone without finding the real optimal solution.The structure of SFS comprises three main stages for producing new solutions, such as the diffusion process, the first update process, and the second update process.In the first stage, each old solution is newly updated more than one time and there are at least two new solutions produced around each old solution.On the contrary, the two remaining stages produce fewer solutions dependent on solution quality in which the worst solutions have more chances to be newly updated while the best solutions have a lower possibility to be newly produced.In MSFS, we concentrate on the improvement of stage 2 and stage 3 by proposing two modifications.In the first modification, two step sizes are suggested to be used and a new technique is proposed for determining one more appropriate step size for updating new solutions.Furthermore, search spaces can be around each individual solution or around the best solution depending on the comparison between fitness function of the considered solution and the best solution.In the second modification, selected solutions to be newly updated are different from those in SFS.In SFS, the worst solutions have more chances to be newly updated rather than the best solutions in the two-update process.The strategy limits to exploit search spaces around solutions with high quality and miss good search spaces.In MSFS, the worst solutions are given priority to be updated in the first update process while the best solutions are selected to be updated in the second update process.In general, MSFS can overcome major shortcomings that SFS has coped with so far, but the implementation of the proposed MSFS method for obtaining promising results has also encountered several difficulties as follows: (i) Set the most appropriate value to the probability of updating bad solutions.The probability is in the range from zero to one, but it should be selected by evaluating the results of the first some cases and then the best value is fixed for other cases.However, a randomly selected value in middle point, such as 0.6 or 0.5, can result in better solutions than SFS method.(ii) Select the best values for population size and iterations.As setting the control parameter to high values, it is sure that the best performance is obtained but execution time may be long.Therefore, the best values must get two advantages, high-quality solutions and short enough simulation time.
The application of these modifications is implemented and evaluated by testing the proposed MSFS method and SFS method on four systems with twelve cases with a different number of generators and different constraints.In addition, the proposed method is also compared to other existing methods.The main contributions of the study can be summarized as follows: (i) Describe the SFS method in detail and point out its advantages (ii) Propose highly effective modifications on SFS in aim to reduce values of population size and iterations, shortening simulation time and finding high performance.(iii) Show the real performance of the proposed method when solving complicated systems.This can help readers to decide if the method should be used for their applied problems.
Other remaining parts of the paper are as follows: Section 2 describes the ELD problem in detail by explaining the objective function and considered constraints.Section 3 introduces the SFS method and then proposes MSFS to the overcome shortcomings of SFS.Section 4 instructs the application of MSFS for solving the ELD problem.Section 5 shows obtained results from MSFS and compares the results with those from other methods including SFS.Section 6 concluded the work in the paper and introduces future work.In addition, Appendix A is also added for showing optimal solutions obtained by MSFS.

Objective Function
The main duty of the ELD problem is to determine the total fuel cost of thermal units and satisfy any constraints of the power system.The total fuel cost can be calculated by using the fuel function of each thermal generating unit that can be expressed as follows:

Traditional Cost Function
When the thermal unit uses single fuel cost, the objective function of ELD can be expressed as a quadratic function and formulated as follows: So, the mathematical modeling for the objective function of the ELD problem is calculated by determining the total fuel cost of N units as follows: where C k is the fuel cost of thermal generating unit k and its form is as shown in Equation ( 2) and in Figure 1.

Objective Function
The main duty of the ELD problem is to determine the total fuel cost of thermal units and satisfy any constraints of the power system.The total fuel cost can be calculated by using the fuel function of each thermal generating unit that can be expressed as follows:

Traditional Cost Function
When the thermal unit uses single fuel cost, the objective function of ELD can be expressed as a quadratic function and formulated as follows: So, the mathematical modeling for the objective function of the ELD problem is calculated by determining the total fuel cost of N units as follows: where Ck is the fuel cost of thermal generating unit k and its form is as shown in Equation ( 2) and in Figure 1.Traditionally, fuel cost function was represented as the second order form once the effects of thermal valves during the power change were not taken into account.However, it is more realistic since the valve effects are considered for the process of increasing or decreasing the power output of thermal generating units.The valve effects are represented as a sinusoidal term.Therefore, the fuel cost function is the sum of the second order Function (1) and the sinusoidal term as shown in Equation (3).The effects can be identified by observing the curve in red in Figure 2 [2].
( ) Traditionally, fuel cost function was represented as the second order form once the effects of thermal valves during the power change were not taken into account.However, it is more realistic since the valve effects are considered for the process of increasing or decreasing the power output of thermal generating units.The valve effects are represented as a sinusoidal term.Therefore, the fuel cost function is the sum of the second order Function (1) and the sinusoidal term as shown in Equation (3).The effects can be identified by observing the curve in red in Figure 2 [2].
Figure 2. Fuel cost function with valve-point effects.

Cost Function with Multiple Fuels and Valve-Point Effects
Usually, thermal generating units can be provided with multiple fuels.So, the cost function can be piecewise curves consisting of two or three second order equations depending on the number of fuels.For this case, Ck can be described in Equation ( 4) and plotted in Figure 3.
Figure 3. Fuel cost function for the case of multi-fuel options.
However, in real operation conditions, the power output of the thermal generating units is controlled by their boilers that are influenced by valve-point effects.Hence, the fuel cost function associated with multiple fuels (MF) options and valve-point effects (VPF) can be mathematically formulated as follows [28]:

Cost Function with Multiple Fuels and Valve-Point Effects
Usually, thermal generating units can be provided with multiple fuels.So, the cost function can be piecewise curves consisting of two or three second order equations depending on the number of fuels.For this case, C k can be described in Equation ( 4) and plotted in Figure 3. Usually, thermal generating units can be provided with multiple fuels.So, the cost function can be piecewise curves consisting of two or three second order equations depending on the number of fuels.For this case, Ck can be described in Equation ( 4) and plotted in Figure 3.
Figure 3. Fuel cost function for the case of multi-fuel options.
However, in real operation conditions, the power output of the thermal generating units is controlled by their boilers that are influenced by valve-point effects.Hence, the fuel cost function associated with multiple fuels (MF) options and valve-point effects (VPF) can be mathematically formulated as follows [28]: However, in real operation conditions, the power output of the thermal generating units is controlled by their boilers that are influenced by valve-point effects.Hence, the fuel cost function associated with multiple fuels (MF) options and valve-point effects (VPF) can be mathematically formulated as follows [28]:

Constraints
Real power balance: The total generated power must be equal to the sum of the total load demand and transmission line losses as the following rule: N k=1 AP k = P LD + P TLL (6) where P LD is the total load demand (in MW).P TLL is total transmission line losses (in MW) and can be calculated by using Kron's formula below: where B kh , B 0k , B 00 are B active power loss matrix coefficients and computed by the following way [44]: (i) Run optimal power flow for a solution to record voltage magnitude and phase angle of all buses (ii) Determine branch currents and impedance matrix (iii) Form Hermitian matrix H (iv) Calculate B coefficients The value of coefficients has a great influence on the total transmission line losses and these coefficients should be updated per iteration due to the change of solutions from optimal power flow program.However, in this study, these coefficients remain unchanged with an assumption that new generation schedule and previous schedule are slightly different or approximately similar [44].
Active power output constraint: Active power output of each generator should be between its upper and lower power limits as shown in Equation (8) below: Prohibited operation zone constraint: In some situations, the thermal units have some prohibited operation zones (POZ) owing to the survival of some weakness in generators or in accessories.These prohibited operation zones create intermittent regions in the cost function.Therefore, when the units operate, it should avoid these regions as shown in the following model:

Stochastic Fractal Search Algorithm (SFS)
SFS was first proposed by Salimi in 2014 for stacking optimization problems [43].SFS algorithm was modified based on basic fractal search (FS) to get rid of the disadvantages of FS, such as: 1.
There are too many parameters which have to be carefully selected and the selection of these parameters must be suitable to avoid trapping in local solutions.

2.
The speed of convergence of the algorithm is too slow because there is no exchange of information between individuals within the group.
SFS can deal with the above drawbacks by proposing two processes: the first update and the second update.So, the structure of SFS consists of the diffusion process, the first update process and the second update process corresponding to three times of newly generated solutions.Concrete information about such SFS is presented as follow:

Diffusion Process
The first process is considered as a narrow zone search technique based on the natural phenomenon of growth.In such process, each solution d from the population (N pop ) will diffuse into the number of diffusions (N d,f ) of new solutions by using Gaussian walk model.This process is carried out by the selection of Equation (10) or (11).
where Gaussian(X best , ∆) and Gaussian(X d , ∆) are two solutions, which are, respectively, found around X best and X d by using Gaussian distribution with the standard deviation of ∆ [43] in which ∆ is obtained by Equation ( 12) below.
where I ter is the current iteration.X best is represented as the position of the best solution and X d is represented as the dth solution in the population.After using both mentioned equations to update new solutions, the quality of all updated solutions X new d, f is measured by determining their fitness function, which is represented as FF d,f (f = 1,..., N d,f ).Then, Fitness function values of new solutions (FF d,f ) are compared with each other to find the lowest fitness and the solution with the lowest fitness, called FF new d and X new d , respectively [45].Finally, the comparison between the newly updated solution' fitness function (FF new d ) and that of the old solution (FF d ) is executed to keep a better one.

The First and Second Update Processes
The first update process is considered a large zone search technique in charge of newly producing solutions around old solutions, which were retained after the diffusion process.First of all, all previous solutions are ranked and assigned to rank d based on their fitness function in which worse solutions are with smaller order and better solutions with higher order.Namely, the best solution with the lowest fitness function is ranked last corresponding to the N pop -th order while the worst one is ranked first.Similarly, other solutions are assigned to their order.In the next step, a ratio (PG d ) is calculated using Equation (13).
In the last step, old solutions receive the decision to update.If the solution owning PG d with a smaller value than a random number θ d , the solution will be newly updated based on Equation (14).We noted that θ d is a random value within the range [0, 1].It is different from a random confidence threshold B in [46].In [46], authors have applied the random confidence threshold B based on different distribution functions, such as beta distribution function, truncated normal distribution function and two-point distribution function to determine the critical confidence threshold E(B) for a continuous opinion model.The value of θ d is used for producing a possibility of updating new solutions rather than producing new solutions with high quality.
In case that PG d is not less than θ d , there is no update action for the considered solution.The whole description can be summarized based on the following Algorithm 1.In the second update process, all steps are similar to those in Algorithm 1 but the Equation ( 14) for updating new solutions is replaced with Equations ( 15) and ( 16) below.

The Newly Proposed Technique for the Two Update Procedures
In the first update procedure of SFS, new solutions via Gaussian walk are updated by using a random solution with step size (∆X) based on a random solution and solution d as Equation (14).Obviously, the search mechanism is a random search way without a good strategy because its generated solutions always search around the random solution with the small step size.So, the first modification of MSFS aims to manage the restrictions of SFS by using a novel search strategy consisting of Equations ( 17) and (18).
Clearly, the search strategies applying two models in Equations ( 17) and ( 18) have a big difference from Equation ( 14) because the proposed model of Equation ( 17) is used to update new solutions around the old solutions whilst the purpose of applying the model in Equation ( 18) is to create new solutions around the so-far best solution.The choice of applying either the model of Equation ( 17) or (18) should be determined by a correct criterion.The criterion is the result of the comparison between EF d and EF aver .The two new definitions are shown in Equations ( 19) and ( 20) as below: where FF d is the fitness function value of each solution X d , FF best is the fitness function value of the current best solution X best .
For the case that EF d is bigger than EF aver , it implies that the current solution is far away from the so-far best solution, new solutions are updated by utilizing Equation (18).On the contrary, it can be understood that the current solution is close to the so-far best solution.Thus, the newly generated solutions are produced by employing Equation (17).
Moreover, the step size has also played a major part in finding a good solution quality.As seen, the step size in Equation ( 14) is always small.Thus, it may be not useful for searching optimal solutions and time consuming for being convergent to the optimal global solutions.To escape the local search zone and handle the mentioned disadvantages, we propound two step size types in order to extend the practicable search zone and overcome the local minima.The details of selecting the step size are presented as follows: where D is the probability of selecting step size types, which can be from 0.1 to 1.By the experience, D should be set by 0.1 to 0.3.This indicates that the step size based on two random points is only used up to thirty percent, while another step size called the large step size can be used up to seventy percent.Consequently, the contribution of the large step size is more significant than that of another.The proposed update process can be easily understood as observing Figure 4.In the Figure, we suppose that there are five solutions in the population.Solution 1 and solution 2 in blue together with the best solution in red are in the group with EF d not higher than EF aver while solution 3 and solution 4 in black are in another group with EF d higher than EF aver .There are two main differences in the process of updating new solutions.Search space around good solution group is exploited while that around bad solution group is neglected.Instead, the best solution is a central solution that is used to find new solutions for bad solutions.In fact, two arrows pointing to the best solution from solutions 3 and 4 illustrates the statement.It is clearly seen that spaces around solution 1, solution 2 and the best solution are marked and sought many times while spaces around solution 3 and solution 4 are abandoned.That is the first difference between good solution group and bad solution group.In the second difference, it can be identified that search spaces around solution 1 are more narrow than those around solution 2 because ε 1 is smaller than random number D but ε 2 is equal to or higher than D. Similarly, space around the best solution is plotted by two circles with different radius, one circle with a higher radius and another with a smaller radius.The smaller circle is the search space for solution 3 because ε 3 is smaller than D whereas the large circle is the search space for solution 4 because D is not higher than ε 4 .The implementation of the proposed update process can be accomplished by using Algorithm 2 below.Algorithm 2. The proposed search strategy for the first update process Calculate given probability of solution d by utilizing Equation (13) if εd < D ( )

The Modifications on Selected Solutions for the Two Update Procedures
Equation ( 13) of SFS indicates that most of the solutions with bad fitness are always updated while solutions with better fitness have a lower possibility to be newly updated.Clearly, PGd of the best solution is the largest value while that of the worst solution is the smallest value.Besides, the value of a random number of solution d, θd is always smaller than one and bigger than zero.Therefore, considering the comparison of PGd and θd, only search spaces around the worst solutions are highly exploited while search spaces around the best solutions are disregarded.The manner has caused the main disadvantage for SFS in seeking solutions in local optimal zones or nearby global optimal solution.Therefore, in the modification, we suggest updating bad quality solutions in the first update procedure while good quality solutions are focused to be newly produced.For implementing the modification, the ratio of updating bad quality solutions Rnbs (where Rnbs is less than one and higher than zero) must be selected in advance and then the number of bad quality solutions is determined and updated in the first update process.Consequently, in the first update, (Rnbs×Npop) worst solutions are newly updated.After the first update, all solutions are ranked and then [(1 − Rnbs) × Npop] best solutions are newly updated in the second update process.
The whole search process of the proposed MSFS method is described in detail in Figure 5 below.

The Modifications on Selected Solutions for the Two Update Procedures
Equation (13) of SFS indicates that most of the solutions with bad fitness are always updated while solutions with better fitness have a lower possibility to be newly updated.Clearly, PG d of the best solution is the largest value while that of the worst solution is the smallest value.Besides, the value of a random number of solution d, θ d is always smaller than one and bigger than zero.Therefore, considering the comparison of PG d and θ d , only search spaces around the worst solutions are highly exploited while search spaces around the best solutions are disregarded.The manner has caused the main disadvantage for SFS in seeking solutions in local optimal zones or nearby global optimal solution.Therefore, in the modification, we suggest updating bad quality solutions in the first update procedure while good quality solutions are focused to be newly produced.For implementing the modification, the ratio of updating bad quality solutions R nbs (where R nbs is less than one and higher than zero) must be selected in advance and then the number of bad quality solutions is determined and updated in the first update process.Consequently, in the first update, (R nbs × N pop ) worst solutions are newly updated.After the first update, all solutions are ranked and then [(1 − R nbs ) × N pop ] best solutions are newly updated in the second update process.
The whole search process of the proposed MSFS method is described in detail in Figure 5 below.

Dealing with the Constraints
In ELD problem, solving load demand-supply balance constraint is one of the most important affairs because it affects the good solution quality and the speed of convergence.So, many methods have used two variable types to deal with the constraints.The first variable type is called a dependent variable and selected from one thermal generating unit.The second variable type is called a control variable and selected from other thermal generating units.In the search strategy of the proposed method, the variables are included in the position of each point in the initialization step and are updated in each iteration.Consequently, the position of point d will include thermal generating unit from unit 2 to unit N as the following formula. where As a result, the balance constraint can be solved by the expressions below.

Handling the Violation of the Dependent Variable
The value of AP 1,d obtained from Equation ( 26) may violate constraint in Equation (8).Therefore, the violation must be evaluated in the quality of solutions.The task is performed by computing the following penalty term.
As seen from Equation ( 27), the thermal generating unit 1 will be penalized if it is higher than upper bound or smaller than lower bound.On the contrary, it will not be penalized or the penalty term will be zero if it is within the upper bound and the lower bound.

Handling the Violation of Upper and Lower Boundaries
After newly updating the power output of from the second unit to the last unit, these units must be checked and corrected exactly by using the following equation:

Handling the Violation of Prohibited Operation Zones
If the active power output of the kth thermal generator falls into one out of prohibited operation zones, it should not be kept but correction and penalty need to be made depending on thermal generators.From the second generator to the last generator (AP k with k = 2, . . ., N) should be corrected by using Equation ( 29) while the first unit AP 1 should not be corrected but it needs to be penalized by using Equation (30).The two equations are as follows: In Equations ( 29) and (30), average powers are defined as the following model:

Handling the Violation of the Spinning Reserve Power (SRP)
Under normal conditions, the operating reserve capacity of the generating thermal units has played an important role in meeting load demands within a short time in case a generator breaks or there is another disruption to the supply.If that is not satisfactory, the approach is penalized.This problem is presented as the following formula:

Calculating Fitness Function
Most of the meta-heuristic methods need to calculate the fitness function of all the solutions in order to arrange the rank of all the solutions.Accordingly, the fitness function is determined as the following equation.

The Detail of MSFS's Procedure for ELD Problem
Step 1. Select parameters for the proposed MSFS, such as the number of points N pop , the maximum value of iterations N Iter , and the ratio of updating bad quality solutions R nbs Step 2. Randomly initialize population based on the rule shown in Equation ( 23 Step 3. Produce new solutions of the diffusion process using Equations ( 10) and (11) Step 4. Correct new solutions by using Equation ( 28) -Handling POZ constraint for unit 2 to unit N based on (29) -Calculate dependent variables AP 1,d by using (26) -Penalize the violation of generation limits for AP 1,d by using Equation (27).
-Penalize the violation of POZ constraint for AP 1,d by using Equation (30).
-Calculate fitness function FF d for the solution d by using (33) Step 5. Compare new solutions and old solutions to retain better ones.
Step 6. Update bad solutions by using the first update process in Algorithm 2 -Correct new solutions by using Equation (28) Step 7. Handling POZ constraint for unit 2 to unit N based on (29) -Calculate dependent variables AP 1,d by using ( 26) -Penalize the violation of generation limits for AP 1,d by using Equation ( 27).
-Penalize the violation of POZ constraint for AP 1,d by using Equation (30).
-Calculate fitness function FF d for the solution d by using (33) Step 8. Compare new solutions and old solutions to retain better ones.
Step 9. Update good solutions by using the second update process in Algorithm 2 -Correct new solutions by using Equation (28) Step 10.Handling POZ constraint for unit 2 to unit N based on (29) -Calculate dependent variables AP 1,d by using (26) -Penalize the violation of generation limits for AP 1,d by using Equation (27).
-Penalize the violation of POZ constraint for AP 1,d by using Equation (30).
-Calculate fitness function FF d for the solution d by using (33) Step 11.Compare new solutions and old solutions to retain better ones.
Step 12. Determine the best solution with the smallest fitness function value.

Numerical Results
In this section, the proposed MSFS and SFS are applied to solve the ELD problem in four different test systems with 3, 6, 10 and 20 thermal units.The first test system considers two cases (cases 1-2), the second system evaluates four cases (cases 3-6), the third system analyzes five cases (cases 7-11), and the final system investigates only one case (case 12).The different cases are due to the different fuel cost function characteristics, different loads and different constraints in which fuel cost functions consider single fuel with and without valve-point effects (VPE) and multi-fuels (MF) with and without valve-point effects while different constraints are power losses (P TLL ), prohibited operating zone (POZ) and power balance.As a result, there are twelve study cases with a detailed description shown in Table 1.
Similar to the meta-heuristic algorithms, the proposed MSFS also needs to set some parameters, such as the population size, the maximum iteration number, the number of diffusions and the ratio of updating bad-quality solutions reported in Table 2.All study cases of the method have been performed in Matlab program language and a computer with 4 GB of Ram and 2.4 Ghz processor.

The Analysis of the Performance Improvement of MSFS
In order to test the level of improvement of MSFS compared to SFS, this section surveys the effect of population and the number of iterations on the results of the methods.The first survey is studied in a complex constraints case, specifically, case 6 with transmission power losses and prohibited operating zones as constraints.The results achieved by SFS and MSFS methods with respect to the best cost, mean cost, worst cost and standard deviation cost are shown in Table 3 and Figure 6.Among the four cost values, the best cost and standard deviation cost are the two most important values considered as major comparison criteria.It means that the comparison of the minimum cost is used to evaluate the best optimal solution and the comparison of the standard deviation cost is used to evaluate the stabilization of finding ability.From Table 3, the population of the two methods is set to five but the number of iterations is changed from 15 to 60 iterations.As observed from the minimum cost, it only needs 20 iterations for MSFS to get the best cost of 15,443.0752($/h).Meanwhile, the best cost of SFS gradually decreases from 15,443.1775($/h) to 15,443.0752($/h) when the number of iterations of SFS is changed from 15 to 60.Specifically, a cost of 15,443.1775($/h) is corresponding to the setting of 15 iterations, a cost of 15,443.1381($/h) is corresponding to the setting of 20 iterations, a cost of 15,443.1197($/h) is corresponding to the setting of 30 iterations, a cost of 15,443.1171($/h) is corresponding to the setting of 45 iterations, and, a cost of 15,443.0752($/h) is corresponding to the setting of 60 iterations.Clearly, to obtain the minimum cost of 15,443.0752($/h), SFS must use 60 iterations but MFSF only uses 20 iterations.As seen Figure 6, the minimum cost of MSFS is 15,443.0772($/h) at 15th iteration, but from the 20th iteration to the 60th iteration, cost of MSFS is also 15,443.0752($/h) and does not change.On the contrary, the best cost of SFS tends to be decreased from 15,443.1755($/h) to 15,443.0752($/h) corresponding to from the 15th iteration to the 60th iteration.This implies that MSF is more efficient than SFS.Besides, MSFS always has a lower standard deviation cost than SFS when the two methods use the same number of iterations.This points out that MSF is more robust than SFS.For other remaining cases, the strong point of the proposed method over SFS is the same manner.The obtained results of MSFS and SFS in regard to the best cost, the average cost, the worst cost and the standard deviation are seen in Table 4 for the rest of the cases.With the same population and iterations, the best cost of MSFS is always smaller than those of SFS.In addition, to see the improvement of MSFS over SFS, we have calculated the reduction cost and added it as the last column in Table 4.For different cases or different test systems, the value of reduction is also different.For example, that of the first test system is from 0.02 ($/h) to 0.04 ($/h), for the second test system is from 0.019 ($/h) to 0.409 ($/h), for the third test system is from 0.0001 ($/h) to 0.0429 ($/h), and for the last test system is 0.006 ($/h).In general, the reduction cost is not high for each hour but the reduction for one day with 24 h or one year with 8760 h is highly significant.Furthermore, the improvement of the proposed method over SFS can be identified more clearly as focusing on average cost, maximum cost as well as standard deviation cost.For other remaining cases, the strong point of the proposed method over SFS is the same manner.The obtained results of MSFS and SFS in regard to the best cost, the average cost, the worst cost and the standard deviation are seen in Table 4 for the rest of the cases.With the same population and iterations, the best cost of MSFS is always smaller than those of SFS.In addition, to see the improvement of MSFS over SFS, we have calculated the reduction cost and added it as the last column in Table 4.For different cases or different test systems, the value of reduction is also different.For example, that of the first test system is from 0.02 ($/h) to 0.04 ($/h), for the second test system is from 0.019 ($/h) to 0.409 ($/h), for the third test system is from 0.0001 ($/h) to 0.0429 ($/h), and for the last test system is 0.006 ($/h).In general, the reduction cost is not high for each hour but the reduction for one day with 24 h or one year with 8760 h is highly significant.Furthermore, the improvement of the proposed method over SFS can be identified more clearly as focusing on average cost, maximum cost as well as standard deviation cost.
To give further evidence of the stability of the proposed algorithm, we also test the distribution of the minimum costs over fifty independent trial runs for case 6 and case 11.Case 6 takes transmission power losses and prohibited operating zones into account and case 11 considers multi-fuels options and valve-point effects as constraints.The best costs of 50 runs obtained by MSFS and SFS are plotted in Figures 7 and 8 for case 6 and case 11, respectively.These figures show a random distribution of 50 values of fitness function values obtained by MSFS and SFS over 50 fruitfully independent runs.It means that after each trial run, there is a red point of MSFS and a blue point of SFS allocated on such curves.These points can be increased or decreased without any rules.The deviation between the peak of SFS and MSFS is very high.The fluctuations of SFS are high since all points of MSFS are approximately lied on a line.Furthermore, it can point out that nearly all points of MSFS have the same cost as the best point but only a few points of SFS are around the best point of MSFS.The manner indicates that SFS hardly ever finds an optimal solution.To give further evidence of the stability of the proposed algorithm, we also test the distribution of the minimum costs over fifty independent trial runs for case 6 and case 11.Case 6 takes transmission power losses and prohibited operating zones into account and case 11 considers multifuels options and valve-point effects as constraints.The best costs of 50 runs obtained by MSFS and SFS are plotted in Figures 7 and 8 for case 6 and case 11, respectively.These figures show a random distribution of 50 values of fitness function values obtained by MSFS and SFS over 50 fruitfully independent runs.It means that after each trial run, there is a red point of MSFS and a blue point of SFS allocated on such curves.These points can be increased or decreased without any rules.The deviation between the peak of SFS and MSFS is very high.The fluctuations of SFS are high since all points of MSFS are approximately lied on a line.Furthermore, it can point out that nearly all points of MSFS have the same cost as the best point but only a few points of SFS are around the best point of MSFS.The manner indicates that SFS hardly ever finds an optimal solution.

Comparison and Discussion
In this part, we will conduct an investigation of the best simulation results of the proposed algorithm and those from other approaches available in the paper.Besides, another comparison formula is concerned to be the number of fitness assessments (NGer), which is obtained by equation ( 34) below: where t stands for the number of generations in every iteration.Depending on the structure of optimization algorithms, they have a difference of the number of generations.For example, CSA and ORCSA technique with two solution generations, t is two while t is one for other technique with one solution generation, such as DE, GA, and PSO methods.For the proposed MSFS, the number of generations in each iteration equals the sum of the number of diffusions and number one ( , 1 df N + ).This formula indicates that an algorithm with lower NGer is assessed to be a more effective method if it has equal or lower best cost.

Test System 1 with 3 Units
In the first test system, a small size system with three generators is dissected with two different study cases in which transmission power losses are investigated in case 1 and valve-point effects are investigated in case 2. The load demand in test system 1 is set to 850 MW.The detailed data input of the first case related to the fuel cost function coefficient and transmission power losses matrix are shown in [17] and [32].In this case, we discuss the results of the proposed method and others such as MOEP [17], NSGA-II [32], ICSA [25], CSA [22], and BBO [15] regarding the best cost, the population (Npop), the number of iterations (NIter), and the number of fitness evaluations (NGer) in Table 5.As seen in Table 5, the solution quality of the MSFS method is the same as that of ICSA [25], CSA [22], and BBO [15] but better than that of MOEP [17] and NSGA-II [32].Besides, NGer of MOEP, NSGA-II and BBO is very high, 150,000 for MOEP, 10,000,000 for NSGA-II and 6000 for BBO while NGer of MSFS is only 135.For that reason, MSFS is capable of searching very good quality solutions for case 1.

Comparison and Discussion
In this part, we will conduct an investigation of the best simulation results of the proposed algorithm and those from other approaches available in the paper.Besides, another comparison formula is concerned to be the number of fitness assessments (N Ger ), which is obtained by Equation ( 34) below: where t stands for the number of generations in every iteration.Depending on the structure of optimization algorithms, they have a difference of the number of generations.For example, CSA and ORCSA technique with two solution generations, t is two while t is one for other technique with one solution generation, such as DE, GA, and PSO methods.For the proposed MSFS, the number of generations in each iteration equals the sum of the number of diffusions and number one (N d, f + 1).This formula indicates that an algorithm with lower N Ger is assessed to be a more effective method if it has equal or lower best cost.

Test System 1 with 3 Units
In the first test system, a small size system with three generators is dissected with two different study cases in which transmission power losses are investigated in case 1 and valve-point effects are investigated in case 2. The load demand in test system 1 is set to 850 MW.The detailed data input of the first case related to the fuel cost function coefficient and transmission power losses matrix are shown in [17] and [32].In this case, we discuss the results of the proposed method and others such as MOEP [17], NSGA-II [32], ICSA [25], CSA [22], and BBO [15] regarding the best cost, the population (N pop ), the number of iterations (N Iter ), and the number of fitness evaluations (N Ger ) in Table 5.As seen in Table 5, the solution quality of the MSFS method is the same as that of ICSA [25], CSA [22], and BBO [15] but better than that of MOEP [17] and NSGA-II [32].Besides, N Ger of MOEP, NSGA-II and BBO is very high, 150,000 for MOEP, 10,000,000 for NSGA-II and 6000 for BBO while N Ger of MSFS is only 135.For that reason, MSFS is capable of searching very good quality solutions for case 1.

Test System 2 with 6 Units
The second test system is established by six thermal units and is divided into four cases from case 3 to case 6.In case 3, case 4, and case 5, the fuel cost function of thermal units is a quadratic function, the data implemented in this situation are obtained from [34] and the required load demand is set to be 800, 1200, and 1800 MW without transmission losses.For case 6, both the prohibited operating zones and transmission losses are considered while the active load demand is set to 1263 MW [33].
Table 7 provides the results yielded by MSFS and other approaches.With regard to the best cost, MSFS gets a cost of 8227.09($/h) for case 3, 11,477.09($/h) for case 4 and 16,579.33($/h) for case 5.These results absolutely suppress FCGA [34] and CGA [34] but share the same position with GA [16], CSA [22], and ICSA [25] for cases 3, 4 and 5. Table 8 layouts the minimum cost obtained by seven approaches.From Table 8, three out of seven methods obtained the same best cost value equal to 15,443.08$.They are KHA [20], CSA [21], and MSFS, respectively.This result is considerably lower than those from BBO [14], NAPSO [33], IPSO [30], and GAAPI [41] by 0.02 ($/h), 0.69 ($/h), 0.92 ($/h), and 6.62 ($/h), respectively.Furthermore, the number of fitness evaluations of MSFS is only 300 for cases 3, 4, 5, and 6 and is smaller than that of the others.Clearly, MSFS is superior to these compared methods.Optimal solutions of the systems are shown in Table A2 in Appendix A.

Test System 3 with 10 Units
In this part, a medium scale power system with 10 generators is employed for five cases from case 7 to case 11.Among these five cases, only case 11 takes valve-point loading effects into account but all the cases use the piecewise quadratic function [11].Four levels of the load demand are examined, which are found in the literature: P DL = 2400 MW (in case 7), P DL = 2500 MW (in case 8), P DL = 2600 MW (in case 9), and P DL = 2700 MW (in case 10).The value of the best cost and the fitness evaluations obtained by the proposed algorithm and ten other methods are listed in Table 9.This table reveals that the best cost of MSFS has equally optimal solution quality with eight compared methods and lower cost than two other ones, such as ELANN [10] and IQPSO [28].However, MSFS is ranked at the first position amongst these methods because of having the lowest value of N Ger with 1.200 while AIS [12] and MPSO [27] have used 3000, DE [5] has used 12.000, IQPSO [28] has used 40.000, and MSFLA, GHS, SFLA-GHS, and SDE [38] have used 9.000.Clearly, MSFS is the standout method with the lowest evaluation for these four cases.Another case (case 11) in such test system is also investigated in order to verify the efficiency of the proposed method where both valve-point loading effects and multi-fuels options are assessed.The complete data are given in [35].The best cost, the population, the number of iterations, and the number of fitness evaluations achieved by eleven methods are arranged in Table 10.As Table 10 shows, the best cost of BBO is only equal to 605.639 ($/h) and is lower than those of other methods.Nevertheless, it is difficult to determine whether the result is accurate.Referring to BBO [14], we see that the best solution of generator 3 is 332.02MW.It is clear that this value is corresponding to the second fuel type but the authors [14] have used fuel cost function coefficients of the third fuel type for calculating the fuel cost of generator 3.This mistake has resulted in a smaller fuel cost than the accurate value.Considering the solution quality aspect, the minimum cost of MSFS obtains 623.834 ($/h), is lower than those from DAA [12], APSO [36], RVM-PSO [35], MSFLA [38], GHS [38], and SFLA-GHS [38] and slightly higher than those from three methods, such as IQPSO [32], DEPSO [36], and SDE [38].However, N ger of MSFS has used 9.000 while IQPSO [32] has used 40.000, and DEPSO [36] has used 25.000.Clearly, N ger of MSFS is many times lower than that of IQPSO [32], DEPSO [36].This shows that MSFS outperforms the other methods for case 11.
Table 10.Comparison of the implemented results for case 7.

Test System 4 with 20 Units
In this paragraph, the proposed method is applied to deal with the ELD for a twenty-unit system with the transmission line losses [9].The performance of MSFS and a number of other approaches is tabulated in Table 11.The most interesting finding in Table 11, a cost of 62,455.616($/h) for HPSO-GSA algorithm is very diminutive but this result is unachievable.Consulting HPSO-GSA [37], we see that the results of the total generated real powers and system power losses do not satisfy the power balance constraints Equation (6).Standing on this view, the cost of MSFS, CSA [25], and ORCSA [25] is only 62,456,633 ($/h) and the value is lower than those of lambda-iteration [9], HM [9], and BBO [31].Moreover, N Ger of MSFS is lower than those of CSA [25], ORCSA [25] and BBO [31].Consequently, MSFS is one of the most effective methods with the lowest fitness evaluation.

Conclusions
In this paper, a revamped version of the conventional stochastic fractal search algorithm, called modified stochastic fractal search algorithm, was proposed for solving the smooth or non-smooth ELD problem with different fuel cost function characteristics, different loads and different constraints.Based on two proposed modifications, the proposed MSFS has become a strong tool with a good solution quality, fast convergence, and stabilization of searching ability.The first modification is in charge of exploiting the global search space while the second one undertakes to enhance solution quality.Therefore, the proposed method has solutions better than SFS.As tested on four systems, fuel cost-saving level over SFS could reach 0.04 ($/h) for a three-unit system, 0.409 ($/h) for a six-unit system, 0.0429 ($/h) for a 10-unit system, and 0.06 ($/h) for a 20-unit system.Furthermore, the performance of the proposed method over other mentioned methods has also been scrutinized.The reduction in cost can be up to 16.65 ($/h) for the second test system, 0.282 ($/h) for the third test system and 0.16 ($/h) for the last test system.In addition, other costs of the proposed method, such as mean cost and maximum cost, are also less than those of the other ones.Consequently, it can be concluded that the proposed method is a useful optimization approach for searching solutions of the ELD problem.
With advantages of the proposed method, in future research, MSFS can be introduced to deal with ELD problems considering complicated models of thermal generating units [47], the combination of heat generators and power generators [48].Besides, the considered ELD problems regarding renewable energies, such as wind power plants and solar power plants, are also interesting studies [49].

Conflicts of Interest:
The authors proclaim that the publication of this paper does not conflict with any personal circumstances or interest.

Figure 1 .
Figure 1.Fuel cost function for the case of single fuel option.

Figure 1 .
Figure 1.Fuel cost function for the case of single fuel option.

Figure 2 .
Figure 2. Fuel cost function with valve-point effects.

Figure 3 .
Figure 3. Fuel cost function for the case of multi-fuel options.

Algorithm 1 .
The first update process of SFS method (i) Sort all solutions by based on the fitness function (ii) Determine rank d value for all sorted solutions (iii) Determine PG d by using Equation (13) For d = 1:N pop if PG d < θ d Update solution d by using Equation (14) else Solution d does not change end end

Figure 4 .
Figure 4.The description of the first update process of the proposed modified stochastic fractal search algorithm (MSFS) method.

Figure 4 .Algorithm 2 .
Figure 4.The description of the first update process of the proposed modified stochastic fractal search algorithm (MSFS) method.

Figure 5 .Figure 5 .
Figure5.The flowchart of using the proposed method for solving a typical optimization problem.
) -Handling POZ constraint for unit 2 to unit N based on (29) -Calculate dependent variables AP 1,d by using (26) -Penalize the violation of generation limits for AP 1,d by using Equation (27) -Penalize the violation of POZ constraint for AP 1,d by using Equation (30).-Calculate fitness function FF d for the solution d by using (33) -Find the best solution X Best with the lowest fitness -Start the first iteration (I ter = 1)

Figure 6 .
Figure 6.Fuel cost obtained by SFS and MSFS with different values of N Iter case 6.

Figure 7 .
Figure 7. Distribution of 50 runs obtained by SFS and MSFS for case 6.Figure 7. Distribution of 50 runs obtained by SFS and MSFS for case 6.

Figure 7 .
Figure 7. Distribution of 50 runs obtained by SFS and MSFS for case 6.Figure 7. Distribution of 50 runs obtained by SFS and MSFS for case 6.

Figure 8 .
Figure 8. Distribution of 50 runs obtained by SFS and MSFS for case 11.

Figure 8 .
Figure 8. Distribution of 50 runs obtained by SFS and MSFS for case 11.

Author
Contributions: L.H.P. and M.Q.D. have coded MSFS for solving systems of ELD problem.L.H.P., M.Q.D. and T.T.N. have written Sections 1, 3 and 4. V.-D.P. and H.-N.N. have been in charge of Section 2, Section 5, Section 6 and other duties.

Funding:
This research received no external funding.

Table 1 .
Summary of twelve study cases.

Table 2 .
The selection of population size and the highest iteration number.

Table 3 .
Results obtained by MSFS and stochastic fractal search (SFS) methods for case 6.
Figure 6.Fuel cost obtained by SFS and MSFS with different values of NIter for case 6.

Table 4 .
Results obtained by MSFS and SFS methods for the rest of the cases.

Table 4 .
Results obtained by MSFS and SFS methods for the rest of the cases.

Table 5 .
Result comparisons for case 1.

Table 6 .
Comparison of the implemented results for case 2.

Table 7 .
Comparison of the implemented results for case 3, case 4, and case 5.

Table 8 .
Comparison of the implemented results for case 6.

Table 11 .
Comparison of the implemented results for case 12.

Table A1 .
Lower and upper limits of the jth prohibited operation zones of the kth thermal generator AP min k , AP max k The upper and lower power output of the kth thermal generator AP min kMk , AP max kMk The upper and lower power output for the fuel type M k of the kth thermal generator AP SRP Sum of spinning reserve power of power system K Amplified factor for the violation.M k Number of the M k fuel type of the kth thermal generator N Penalty term for the violation of thermal generator 1. PEN POZ,AP 1 Penalty term for violating POZ constraint of the first thermal generator PEN SRP Penalty term for the violation of spinning reserve power.vk, x k , y k , s k , t k Constant fuel cost function coefficients of the kth thermal generator v kM k , x kM k , y kM k ,s kM k , t kM k Constant fuel cost function coefficients for the M k fuel type of the kth thermal generator X new1Newly updated position of X d X t1 , X t2 , X t3 , X t4Randomly picked solutions from current population α, β, γ, ε, λ, D Random number with value in the range [0, 1] The best solutions for system 1 found by MSFS. d

Table A2 .
The best solutions for system 2 found by MSFS.

Table A3 .
The best solutions for system 3 found by MSFS.

Table A4 .
The best solution for system 4 found by MSFS.