Maximizing Total Profit of Thermal Generation Units in Competitive Electric Market by Using a Proposed Particle Swarm Optimization

: In the paper, a proposed particle swarm optimization (PPSO) is implemented for dealing with an economic load dispatch (ELD) problem considering the competitive electric market. The main task of the problem is to determine optimal power generation and optimal reserve generation of available thermal generation units so that total profit of all the units is maximized. In addition, constraints, such as generation limit and reserve limit of each unit, power demand and reserve demand, must be exactly satisfied. PPSO is an improved version of conventional particle swarm optimization (PSO) by combining pseudo gradient method, constriction factor and a newly proposed position update method. On the other hand, in order to support PPSO to reach good results for the considered problem, a new constraint handling method (NCHM) is also proposed for determining maximum reserve generation and correcting reserve generation. Three test systems with 3, 10 and 20 units are employed to evaluate the real performance of PPSO. In addition to the comparisons with previous methods, salp swarm optimization (SSA), modified differential evolution (MDE) and eight other PSO methods are also implemented for comparisons. Through the result comparisons, two main contributions of the study are as follows: (1) NCHM is very effective for PSO methods to reach a high success rate and higher solution quality, (2) PPSO is more effective than other methods. Consequently, NCHM and PPSO are the useful combination for the considered problem.


Introduction
Economic load dispatch (ELD) is one of the most important problems in power systems due to its significant contributions to economy and operation stabilization of power system. The ELD problem is mathematically formulated by the objective of minimizing fuel cost and a set of constraints regarding thermal generation units and power systems [1]. The earliest ELD problem did not consider power losses in transmission lines due to effects of resistance and reactance of transmission lines and also ignored valve effects on power increase and decrease process of thermal generation further investigation of performance, researchers can select the most appropriate option for their own problem in the power system or other engineering fields. In this paper, we apply PSO variants consisting of conventional PSO [33], CF-PSO [39], IW-PSO [40], PG-PSO [41], IW-PG-PSO, CF-PG-PSO, TVIW-PSO [42], TVAC-PSO [43], and PPSO. In PPSO method, we have combined pseudo gradient method, constriction factor, and a newly proposed velocity update method. Pseudo gradient method is useful in determining better direction for moving to new positions meanwhile constriction factor can support to limit search space. The two applications can form a high-performance method, which is CF-PG-PSO. Furthermore, we also propose a new position update method by using the sofar best position of each particle instead of using the previous position in PSO. In addition, we propose NCHM that can satisfy power reserve demand easily and successfully but the effectiveness is high. For presenting the real performance of NCHM, the survey of obtained results from PPSO and eight other PSO methods with and without using NCHM is accomplished by using a three-unit system with convex fuel cost function, a ten-unit system with convex fuel cost function and a twentyunit system with nonconvex fuel cost function. For indicating the real performance of PPSO, it is compared to these eight PSO methods, SSA, MDE and other previous methods such as PSO, DE, CSA, ELF-HNM and five LF-HNN methods. In summary, the novelties of the study are as follows: (1) Propose the effective NCHM for handling constraints.
(3) Consider valve effects on thermal generation units for the considered problem.
In addition, the application of PPSO method also has some advantages as follows: (1) PPSO method has few control parameters, population and the number of iterations. Therefore, the setting of the two parameters is simple. (2) The process of evaluating solution quality is easily and simply performed by calculating fitness function.
By owning novelties and advantages above, the study can reach the following main contributions: (1) Reach very high success rate with 100%: Implemented methods using NCHM always reaches all successful runs but the same implemented methods without using NCHM must suffer much lower than 100% for success rate. (2) Converge to high quality solutions: NCHM supports implemented methods to find global optimum solutions with fast speed and reach high stability. (3) PPSO method always reaches better results than other PSO, SSA, MDE and previous methods. (4) PPSO method is faster than approximately all other methods for study cases.
However, in order to reach good results and appropriate simulation time for PPSO method, some difficulties exist as follows: (1) The most appropriate values for the population and the number of iterations are not easy to select. In fact, higher values can result in better results but simulation time is still increased correspondingly. If high values are set, all methods have the same best solution and the evaluation is not exactly performed. In this case, real performance of PPSO method cannot be shown. (2) The procedure of applying PPSO method is a long iterative algorithm. Therefore, the implementation procedure must be careful and verification procedure must be serious.
In addition to introduction, other main sections of the paper are as follows: Mathematical formulation of the problem is shown in Section 2. The construction of applied PSO methods is expressed in Section 3. The application of the PSO methods for the problem is described in detail in Section 4. Comparison and discussion of obtained results are shown in Section 5. Finally, a summary of the contribution and future work is included in Section 6.

Objective Function
The ELD problem in the competitive electric market is established by the presence of an objective function and a set of constraints regarding thermal generating units as well as power systems. In order to present the considered objective function, the fuel cost function for generating electricity is first mentioned as follows: For the case without considering valve point loading effects on thermal generation units during the electricity generation process, fuel cost function is approximately expressed as the second order function below [12]: However, it is more appropriate since the valve effects are considered for the operation process of the units, and a more complex function is used as follows [28]: In the competitive electricity market, reserve power is really necessary and thermal units must generate power higher than predetermined demand if reserve is allocated and used by loads. So, fuel cost for both reserve generation and power generation must be taken into account based on Formulas (1) and (2) above. The two new fuel cost functions are as follows: sin n n n n n n n n n n n n n In the competitive electric market, each thermal power plant concerns two cases, the payment for power delivered and the payment for reserve power. For each case, total revenue and total cost are calculated as follows: (1) Payment for power delivered In this case, reserve power is only paid if the reserve power is used by customers. Thus, the price of reserve power (PriceRP) is higher than the price of delivered power (PriceDP). The total revenue (TR) and total cost (TC) are calculated by [19]: In the two equations above, r is the probability that reserve power is called by the requirement from loads. Stn is on/off operation status of the nth thermal generation unit. The parameter has two values only, 1 for on operation status and 0 for off operation status. In addition, FSUn is start-up fuel cost of the nth thermal generation unit. The determination of the status and the start-up fuel cost are really important for finding optimal generation and optimal reserve as unit commitment problem is considered. However, in the study, we consider pure economic load dispatch problem and the assumption is that all thermal generation units are working. Therefore, the start-up fuel cost and on/off operation status can be neglected, and the result does not influence the task of determining optimal generation and optimal reserve.
(2) Payment for reserve power In this case, thermal power plants receive the price of reserve power because the reserve power is not used. As the reserve power is used by loads, thermal power plants receive the price for delivered power that was generated. So, the price of reserve power is much less than that of delivered power. TR and TC are obtained by [19]:   In Equations (7) and (8), PriceDP and PriceRP are not fixed values and they are different for different time periods in a day [44,18,19]. In fact, this issue was demonstrated [18] and different values for these prices were then applied [19]. In this paper, we consider only one period for pure economic load dispatch problem. Therefore, we only considered one fixed value for PriceDP and one fixed value for PriceRP for each study case.
Finally, the total profit (TP) is determined by [18]: The objective function is to maximize the total profit or minimize the minus total profit as follows: In the Equations (5) and (7), total reserve power from all units ( 1 N n n RG   ) can be from zero to reserve demand (RD) depending on the obtained total profit. As seen from Equation (10), total profit can be higher than zero if the total revenue is higher than total fuel cost. It can be seen from the total fuel cost functions (6) and (8), fuel cost of generation and reserve can be considered if reserve is called and used. So, Equations (5) and (7) can be obtained for any values of reserve and the profit is dependent on the revenue for selling power and total fuel cost for generating power.

The consisdered Constraints
(1) Power demand Basically, the ELD problem is about determining power generation of all thermal units for minimizing total cost of these units. Power balance constraint in the problem is about active power in which total generation, power loss and power demand must follow the equality constraint [45,46]. Frequency stability is really important in high voltage networks [47] because power system cannot be stably working if the frequency oscillation happens. Thus, in traditional economic load dispatch problem, active power balance is seriously constrained as the following model [45]: where Ploss is the total power losses in all transmission lines; TPG is total power generation of all power sources [46] that can supply electricity to loads. TPG can be expressed in detail as follows: where PHi is the power generation of the ith hydropower plant; PPj is the jth power generation of photovoltaic system; PWk is the power generation of the kth wind turbine; and PGn is the power generation of the nth thermal generation unit; I is the number of hydropower plants; J is the number of photovoltaic systems; and K is the number of wind turbines However, as considering competitive electricity market for economic load dispatch problem, the power balance constraint above is rewritten as follows [18,19]: On the contrary to the traditional economic load dispatch problem (neglecting competitive electric market), total power generation from N units in the problem can be lower than the forecasted power demand as long as the total profit of N units is high. Generation companies can select to purchase power lower than demand [19].
However, power generation of each thermal generating unit must follow the limits below [30] n n n LB PG UB   (14) (2) Reserve power demand In the current industry, the significance of the total mean electricity cost can be directly found in the requirement to satisfy the (N-1) contingency [48]. However, two high changes in the competitive industry that are taking place can lead to the high deviations from the operation concept [18]. In the studies [32,49], the reasons for making reserve markets and reserve generation contract were presented. First, power energy providers will sell power energy to their customers via signed contracts. In order to avoid the service interruption and compensate customers for the interruption due to the failure of their own generator, the power energy providers may buy reserve power energy from other providers. Second, since industrial zones or other high-power energy customers usually move to deregulation of use of electricity, reserve market is really necessary. Namely, the reserve demand constraint is formulated by:  (15) Similar to power generation, total reserve power of all units can be lower than reserve demand as long as total profit is high and generation companies can purchase reserve lower than the forecasted reserve demand [19].
However, reserve capability of each thermal generating unit is not infinite and constrained by [18,19] As showing constraints (16) and (17) above, reserve generation limits of each thermal generating unit are not fixed for different reserve demand and power demand. The limits are dependent on power generation, which was predetermined. This issue is totally different from power generation limits that are shown in constraint (14). So, the study focuses on a good method to reach the most appropriate limits for reserve generation and find the best reserve power for each unit in competitive electricity market. In general, previous studies have ignored the major issue.

CF-PSO and IW-PSO
In 1995, Kennedy and Eberhart [33] developed conventional PSO to solve the optimization problem and then the method was modified to be more effective for more complicated problems [39][40][41][42][43]. PSO method is represented as two typical terms including velocity and position. The two main factors are formulated by: The process of updating new velocity by using (18) was considered to be limited due to the constant change of velocity [39,40]. Then, inertia weight factor [39] and constriction factor [40] were proposed for finding out more promising velocity. The application of the two factors results in the following model: where: In CF-PSO and IW-PSO, position update still utilizes Formula (19).

TVIW-PSO and TVAC-PSO
The combination of inertia weight factor and constriction factor has been applied to develop TVIW-PSO method [42]. The proposal developed the new velocity as follows: On the contrary, the study [43] has proposed the change for two acceleration coefficients c1 and c2, and the new velocity of TVAC-PSO is determined by: In the two equations above, c1initial and c1End are initial and final cognitive acceleration factors, respectively. c2initial and c2End are initial and final social acceleration factors, respectively. The factors are predetermined and fixed during the application of TVAC-PSO for a typical optimization problem.

PG-PSO
The pseudo-gradient based search algorithm (PGBSA) [41] has been proposed with intent to select a more appropriate velocity direction for each considered particle. PGBSA could enable PSO to determine the best direction in a large search space without requesting any complicated computation process. The application of PGBSA for PSO has developed PG-PSO and it has been implemented for nonconvex economic load dispatch problem with complicated constraints and non-differentiable objective function [41]. The method for updating new velocity in PG-PSO still applied formula of PSO meanwhile the method for updating new position is determined by:

The Proposed PSO Method
In this paper, we suggest applying the combination of PG-PSO with inertia weight factor to form IW-PG-PSO and the combination of PG-PSO with constriction factor to form FC-PG-PSO. In addition, we suggest one more modification for updating new positions as follows: By applying the combination of FC-PG-PSO and the newly updated position model, PPSO is first introduced in the paper. In summary, PPSO method applies Formulas (21) and (23) to update new velocity and then applies Formulas (29) and (30) to update new position.

The New Constriaint Handling Method for Reseve Power
As shown in Equation (14), reserve power of each thermal generating unit n can be from 0 MW to (UBn-LBn) but it is totally different if Equation (17) is considered as another main constraint. In fact, Equation (17) indicates that the reserve must not be higher than (UBn-PGn) while PGn can be higher than LBn. So, we suggest that upper bound of reserve should be determined first and then real reserve can be constrained by the upper bound. Namely, the two formulas below should be used.
Therefore, in order to sure reserve power always satisfy constraints (16) and (17), the upper bound of reserve should be determined first by using Equation (31) and then its value is checked and corrected by the following model:

Selection of Control Variables and Population Initialization
Basically, each solution p contains a set of control variables corresponding to considered problems. In the paper, power generation and reserve generation are selected to be control variables and included in each solution p (Pop) as the following expression:

Calculation of Dependent Variables
After having power generation and reserve generation from unit 2 to unit N, the power generation and the reserve generation of the first thermal generating unit is determined by

Correction for Produced Control Variables
After initializing solutions by using Equation (37), reserve generation must be checked and corrected if it is outside the range between lower bound and upper bound. The correction for reserve can be accomplished by using Section 4.1.
In addition, after updating new solutions by using Equation (30), power generation must be also corrected by using Equation (40) Finally, reserve generation continues to be corrected by using Section 4.1.

Handling Violation of Power Demand and Reserve Demand
After determining power generation and reserve generation by using Section 4.2.2, power demand and reserve demand are checked and penalized if violations happen. Two penalty terms corresponding to the violation of power demand and reserve demand are calculated by using the following models:

Handling Violation of the First Thermal Generating Unit
Power generation and reserve generation of the first thermal generating unit cannot be corrected because they have a huge contribution to exactly met power demand and reserve demand. Thus, penalty terms for the violations of the first thermal generating unit must be taken into account and determined by:

Fitness Function
Fitness function must be calculated to evaluate the quality of solutions. Hence, the fitness function has to point out the quality of objective function and the violation level of constraints and dependent variables. In an ELD problem with the competitive electric market, the fitness function is the sum of minus total profit and penalty terms for the violations of power demand, reserve demand and the first thermal generation unit. In case of finding a valid solution, fitness function and the minus total profit are the same meanwhile penalty terms are zero. Nevertheless, there is no warrantee that the valid solution is the global optimum or close to the global optimum. In the study, the fitness function is expressed as follows: where K1, K2, K3 and K4 are penalty factors and determined by experiment.

Establishing Limits of Velocity and Producing Initial Velocity
The upper limit and lower limit of velocity are respectively determined by: where Lim is the velocity limit factor and can be selected from 15% to 20% [36]. Similar to initial solutions, initial velocity has to be determined by:

Termination Criterion for Iterative Algorithm
Generally, termination criteria for an iterative algorithm in solving optimization problems can be maximum mismatch of considered constraints, error tolerance of two consecutive iterations or the number of iterations dependent on characteristic of implemented methods and characteristic of considered problems. In the study, PPSO is a population-based method and heavily influenced by random factors. Consequently, the number of iterations is used for the stopping criterion and collecting results. A high possibility is that the high number of iterations can result in better optimal solutions and more stable search ability. Thus, the search process is carried out until the last iteration is reached. The most appropriate selection of the iteration is reached by experiment. However, the performance of a run is heavily influenced by iterations and population. PPSO handles constraints (13) and (15) by selecting reasonable control variables shown in (34) and other constraints (14), (16)-(17) by correcting and using other penalty methods in Equations (31)- (44) so that considered constraints are exactly met. One solution is considered as a valid one if there is no violation of the considered constraints. Normally, one solution found by using the optimization algorithm, which is called the optimal solution, does not necessarily have good quality. For an implemented run, if the total profit is equal to fitness function at the last iteration (i.e., G = Gmax), it is a successful run. On the contrary, if fitness function is higher than total profit (corresponding to penalty terms in (41)-(44) are higher than zero), there is at least one violated constraint and this is an unsuccessful run.

The Entire Search Process of PPSO for the Considered Problem
The whole search process of PPSO can be summarized in Figure 1 and described in the following computation steps.
Step 1: Set value to Nop and Gmax for the proposed method.
Step 6: Determine penalty terms for the violation of power demand, reserve demand, power generation of the first thermal generating unit and reserve generation of the first generating unit by using Equations (41)-(44).
Step 7: Determine fitness function for each solution by using (45).
Step 8: Among current solutions, determine the best solution with the lowest fitness function and set to PoGbest.
Step 9: Set current solutions to Pobest,p and current iteration G to 1.
Step 15: Determine penalty terms for the violation of power demand, reserve demand, power generation of the first thermal generating unit and reserve generation of the first generating unit by using (41)-(44).
Step 16: Determine fitness function for each solution by using (45).
Step 17: Compare Pobest,p and new p Po (p = 1, …, Nop) to keep better one, and set kept one to Pobest,p.
Step 18: Among current solutions Pobest,p, determine the best solution with the lowest and set to PoGbest.
Step 19: If G = Gmax, stop search process. Otherwise, set G = G + 1 and back to Step 10.  Figure 1. The flowchart of solving the considered problem by using proposed particle swarm optimization (PPSO) method.

Numerical Results
In this section, we have implemented nine methods including PSO, CF-PSO, IW-PSO, PG-PSO, IW-PG-PSO, CF-PG-PSO, TVIW-PSO, TVAC-PSO and PPSO for three test systems with two different cases. The summary of the test systems and the two considered cases is as follows: Test system 1: Three units with convex fuel cost function shown in Equation (1) Test system 2: Ten units with convex fuel cost function shown in Equation (1) Test system 3: Twenty units with nonconvex fuel cost shown in Equation (2) Case 1: Total revenue and total fuel cost are obtained by using Equations (5) and (6) Case 2: Total revenue and total fuel cost are obtained by using Equations (7) and (8) The whole data of the tests and parameters corresponding to the two cases are given in Tables A1-A4 in the Appendix. All implemented methods are coded on Matlab program language-version R2016a and run on a personal computer with configuration as follows: CPU: Intel Core i7 with 2.4 GHz processor and 4 GB of RAM 4 GB, GPU: Intel HD Graphics 5500, and system version: Windows 8.1 Pro-64-bit. For each study case, 50 successful runs are obtained. Basic parameters of PSO methods are selected as follows: (1) Nop = 5 and Gmax = 5 for test system 1 (2) Nop = 20 and Gmax = 100 for test system 2 (3) Nop = 30 and Gmax = 500 for test systems 3

The Impact of the Proposed NCHM on Results
In this section, we have run nine methods with and without using NCHM. Obtained results that are used for comparison are success rate (SR) of reaching 50 successful runs, the maximum total profit (MTP) and the average total profit (ATP). MTP is the best profit over 50 successful runs meanwhile ATP is the average profit of 50 successful runs. MTP is used to evaluate the ability of finding the best optimal solution meanwhile ATP is used to reflect the stability of the method over 50 successful runs. SR is compared to reflect the ability of dealing with all constraints of applied methods. Normally, methods with higher MTP are more effective because it can find better solutions; however, real improvement of the methods can be further investigated as considering ATP for comparison. In fact, ATP is the average value of 50 successful runs and higher ATP is corresponding to higher quality of 50 successful runs. ATP is more valuable for a small-scale system like the three-unit system since MTP of methods with and without CHM is not highly different, it is even considered approximately the same. However, ATP of methods with and without CHM is much different. In this case, methods with higher ATP is more effective and more stable for finding optimal solutions over a number of successful runs. So, ATP is really necessary for comparison and evaluation. Furthermore, for a more exact comparison, we also calculate higher MTP and higher ATP that methods using NCHM can reach as compared to the same methods without using NCHM. Then, the values are converted to percent, which is similar to the improvement level of NCHM.
As a result, all methods with using NCHM can reach 100% for success rate, but the success rate from these methods without using NCHM is only from 83.3% to 92.5%. Figures 2 and 3 show MTP and ATP for case 1 and case 2 of system 1. Similarly, Figures 4 and 5 show MTP and ATP for the two cases of system 2, and Figures 6 and 7 show MTP and ATP of the two cases of system 3. In these figures, the bars in blue and red are MTP values with and without using NCHM, while the bars in grey and yellow are ATP values with and without using NCHM. It is seen that MTP with and without using NCHM is approximately equal for system 1, while MTP with using NCHM is higher than MTP without using NCHM for system 2 and system 3. On the contrary, ATP with using NCHM is always much higher than ATP without using NCHM for two cases of three systems. The comment is the same for MTP and ATP of nine implemented methods.
Tables 1-3 report better MTP and ATP (in $/h and in %) of methods using NCHM as compared to the same methods without NCHM. Table 1 sees that using NCHM can reach higher MTP from $0.015 to $1.450 for case 1, and from $69.056 to $179.722 for case 2 of system 1. The higher value is similar to the improvement level from 0.001% to 0.132% for case 1 and from 0.037% to 0.133% for case 2. It is clear the MTP is not much improved by using NCHM, but the improvement of ATP is much more significant. In fact, better value of ATP is from $69.056 to $179.722 for case 1 and from $142.252 to $314.008 for case 2. The values are equivalent to the high improvement from 7.347% to 22.305% for case 1 and from 15.825% to 42.9% for case 2. Table 2 shows that MTP can be higher from $39.15 to $938.57 corresponding to the improvement level from 0.27% to 6.89%, and ATP can be higher by from $1741.05 to $4117.138 corresponding to the improvement level from 13.98% to 41.75% for case 1 of system 2. For case 2, MTP can be higher by from $109.84 to $531.67 corresponding to the improvement level from 0.81% to 4.06%, and ATP can be higher by from $9357.84 to $13,026.5 corresponding to the improvement level from 1544.92% to 11,313.34%. Table 3 shows that MTP can be higher up to $425.54 corresponding to the improvement level of 2.09% and ATP can be higher up to $2263.72 corresponding to the improvement level of 12.55% for case 1. For case 2, MTP can be higher up to $286.23 corresponding to the improvement level of 1.97% and ATP can be higher up to $2376.208 corresponding to the improvement level of 20.24%.
In summary, the application of NCHM can support methods to reach significant achievements as follows: (1) Methods using NCHM can reach the highest SR with 100% but SR of the methods without using NCHM is much lower, only from 83.3% to 92.5%. (2) NCHM can support methods to find the global optimum solutions, high search stability and low possibility to low quality solutions.

Comparison for Test System 1
In this section, we compare the real performance of the proposed PSO method with other PSO and previous methods as testing on three-unit system. In addition, we have also implemented salp swarm algorithm (SSA) [50] and modified differential evolution (MDE) [51]  Clearly, PPSO can find the best optimal solution with the highest profit and all runs of PPSO reach higher performance than other PSO methods, SSA and MDE. Furthermore, PPSO is a more stable method in the searching process since its standard deviation (STD) is much lower than these methods. Namely, it is 96.4 for case 1 and 97.3 for case 2 whereas that of other ones is increased from 101.6 to 196.7 for case 1 and from 110.1 and 212.5 for case 1.
As compared to other remaining methods such as ELF-HNM [30], PSO [31], CSA [31], DE [31] and five Hopfield Lagrange network-based methods [31], PPSO can reach the same best solutions as approximately all these methods excluding LF-HLN-GdF [31], LF-HLN-GF [31] and LF-HLN-LF [31] for case 2 with worse solutions than PPSO. In connection with the comparison of ATP, PPSO is more effective than other metaheuristic algorithms like PSO and DE but the achievement is not reached again as comparing to CSA and five Hopfield Lagrange network-based methods. However, PPSO is still superior to these methods since it has been run by setting 5 to population and 5 to iterations and spent less computation time. In fact, DE, PSO and CSA have been implemented by using 5 for population and 500 for the number of iterations meanwhile these five Hopfield Lagrange networkbased methods are deterministic algorithms with very small change among different runs. Figures 8 and 9 show the best convergence characteristics corresponding to the best run over 50 successful runs while Figures 10 and 11 illustrate the mean solution searching characteristic of the 50 successful runs for case 1 and case 2, respectively. The four figures have the same manner that PPSO is slower than other ones at the first iterations but PPSO converges to better solutions at final iterations. Figures 8 and 9 are about the maximum total profit, so all methods have the same final point. On the contrary, Figures 10 and 11 are about the average total profit, so the final point of PPSO is much higher than that of other ones. It is clear that PPSO has better ability of jumping out local zones and converge to more promising zones.
In summary, PPSO can improve result more effectively than other PSO methods and it can reach approximately equal or better results than other ones, however, it is always faster than other ones. Therefore, PPSO is really a highly efficient method for the system with 3 units and convex fuel cost function.  [31] 1102       Figure 11. Mean solution searching characteristic of implemented methods for case 2 of system 1.

Comparison for Test System 2
In this section, we compare the real performance of the proposed PSO method with other PSO, SSA and MDE and previous methods by using the ten-unit system without valve effects of thermal units. Tables 6 and 7 show result comparison for case 1 and case 2 of the system. From the tables, it is seen that PPSO can reach better MTP and ATP than other implemented PSO, SSA and MDE methods for the two cases. MTP and ATP from PPSO are respectively $14,564.74 and $14,193.08 for case 1 and $13,635.12 and $13,525.28 for case 2 whereas those from other implemented methods are much worse. In fact, SSA must suffer the worst MTP with the lowest values of $14,370.95 for case 1 and $13,597.06 for case 2 whereas the second-best methods consisting of CF-PSO and TVIW-PSO can reach $14,563.77 for case 1 and the second-best method, IW-PG-PSO, can reach $13,635.04 for case 2. Similarly, ATP of the worst method and the second-best method is $13,918.55 and $14,128.56 for case 1 and is $13,086.99 and $13,454.77 for case 2. Clearly, all methods cannot reach the highest profit and find the global optimal solution that PPSO can. Furthermore, the stability of PPSO is always better since the standard deviation is also much lower than other ones. The standard deviation of PPSO is 236.9 for case 1 and 105.1 for case 2 whereas that of others is from 237.1 to 617.6 for case 1 and from 109.8 to 605.7 for case 2. The convergence characteristic for the best run and mean solution searching characteristic for 50 successful runs are respectively plotted in Figures 12-15 for the two cases. As seen in Figure 12 and 13, PPSO cannot reach better solutions than other ones at the first 20 iterations, however, PPSO can find better solutions. Moreover, the superiority of PPSO over other ones can be clearly seen through Figures 14 and 15. Mean convergence curves of PPSO have much higher profit than those of other ones for both case 1 and case 2 from the first iteration to the last one. This means that PPSO has stronger search ability than other ones.
As compared to other remaining methods [30,31], PPSO can reach better MTP than approximately all methods for the two cases excluding LF-HLN-EF, LF-HLN-THF and ELF-HNM for case 1. Especially, as comparing to PSO [31], DE [31] and CSA [31], PPSO can reach much higher MTP and ATP. PPSO can reach higher MTP than PSO, CSA and DE by $382.55, $0.69 and $511.71 for case 1, and $1406.67, $929.64 and $1471.55 for case 2. Similarly, PPSO can also reach much higher ATP than PSO, DE and CSA for the two cases. Clearly, the improvement of PPSO over PSO, DE and CSA is significant. Although PPSO can reach higher ATP than all methods, PPSO has implemented only 100 iterations with population of 20, whereas, PSO, DE and CSA have used 500 iterations and population of 5. Furthermore, computation time of the proposed method is still much faster than these methods.
In summary, PPSO can find better or the same solutions with other compared methods but it outperforms these methods in terms of convergence speed and stability of searching ability. Consequently, PPSO is really effective for the system with ten units and without valve effects of thermal units.   Table 7. Comparison of results obtained for case 2 of system 2.

Comparison for Test System 3
In this section, PPSO is compared to eight other PSO methods, SSA and MDE by employing a 20-unit system with valve effects on thermal generation units. MTP and ATP obtained by all implemented methods are respectively plotted in Figures 16 and 17 Figure 20 showing standard deviation of 50 successful runs is also a good evidence for confirming the strong search of PPSO since two bars of PPSO for the two cases are the lowest among eleven implemented ones.
As seen from the computation time shown in Figure 21, all applied methods have approximately equal time because we have set the same population and the same number of iterations for them. So, simulation time is approximately the same. However, the convergence characteristics of the best run and the mean solution searching characteristics in Figures 22-25 for case 1 and case 2 can indicate that PPSO is faster and more stable than other ones. PPSO can find much better solutions than others after the 100th iteration and even solution of other methods at the 500th iteration is much worse than that of PPSO at the 100th iteration. Clearly, PPSO is significantly faster than these compared methods.
In summary, as implementing PPSO and other meta-heuristic algorithms for the largest system with 20 units considering valve effects, PPSO can show outstanding performance, since it can reach much better solutions, more stable searching ability and faster search process. Hence, PPSO is a promising method for dealing with ELD problem considering competitive electric market and valve effects on thermal generation units.         Figure 23. Convergence characteristic of implemented methods corresponding to the best run for case 2 of system 3.  Power generation and reserve power of each thermal generation unit are reported in Table A5, Tables A6 and A7 in Appendix for the three studied systems. As calculated from the tables, total power generation and reserve are, respectively, 924.5042 MW and 100 MW for the three-unit system,  2) for the twenty-unit system. As compared to power demand and reserve demand shown in Table A4, balance of generation and reserve for ten-unit system is satisfied, while the balance of the three-unit system and the twenty-unit system is not met. The same results can be seen as referring to three-unit system and twenty-unit system of unit commitment problem [19], and three-unit system of economic load dispatch problem [30]. In competitive electric market, generation companies can provide power and reserve lower than forecasted power and reserve demand as long as they reach high profit [19].

Conclusions and Future Work
In the paper, a proposed particle swarm optimization has been compared to SSA, MDE and eight other PSO methods in finding optimal solutions of ELD problem taking into account competitive electric market. Study cases were three different systems with 3, 10 and 20 units in which the 20-unit system has considered valve point effects on thermal generation units. In addition, a new constraint handling method has been also applied for all methods. As a result, the proposed constraint handling method was very useful in reaching the highest success rate of 100% and finding much better optimal solutions for all cases. As compared with SSA, MDE and eight other PSO methods, the proposed method was the best because it could find the same or better solutions than these methods but it was faster than the methods for approximately all study cases. Furthermore, the proposed method was also compared to other previous methods for evaluating clear improvement level. The proposed method could find either equal or better solution quality than others meanwhile the proposed method has used smaller number of iterations. As a result, it is recommended that the new constraint handling method should be used for the problem as applying metaheuristic methods and the proposed method should be used for ELD problem considering competitive electric market.
In this paper, we have considered only thermal power plants in competitive electric market, namely thermal generation units in a thermal power plant. It is obvious that all types of power plants can supply electricity in a competitive electric market; however, their optimization operation strategies can be different from thermal power plants. In fact, hydropower plants can store water by using pumped storage system and reaching the best operation strategy. Wind turbine can adjust power output by changing bitch angle. Operating power plants in competitive electric market is separated and different power plants have different prices depending on type of power plants and different fuel characteristics. In the future, we will consider only wind turbines or photovoltaic systems in the electric market. The proposed PSO method will be successfully applied for the renewable energies in electric market. As shown in Section 4, the proposed PSO method can be successfully implemented for the problem as long as fitness function is correctly established. In the fitness function of the considered problem, objective function is total profit meanwhile penalty terms are to avoid the violations of power generation and reserve power of the first thermal generation unit, the violation of power demand and the violation of reserve demand. For the case that wind turbines together with photovoltaic systems are considered and mathematical formulation is successfully developed, the proposed PSO and other meta-heuristic algorithms are capable of solving the new problem. Funding: This study is funded through the Project of the Year 2020 of the Ministry of Education and Training -Vietnam. We also thank Ho Chi Minh City University of Technology and Education for the support.

Conflicts of Interest:
The authors declare that there is no conflict of interests regarding the publication of this paper. Abbreviation IW-PG-PSO Inertia weight factor and pseudo gradient -based particle swarm optimization CF-PG-PSO Constriction factor and Pseudo gradient-based particle swarm optimization TVIW-PSO Time varying inertia weight factor-based particle swarm optimization LF- Power generation of the first thermal generation unit Pobest,p The so-far best position of the pth particle PoGbest The so-far best position of all particles RD Forecasted reserve power demand RGn Reserve generation of the nth thermal generating unit Old velocity and position of the pth particle ω Inertia weigh factor ωmin, ωmax Minimum and maximum value of inertia weigh factor