Predictor-Corrector Primal-Dual Interior Point Method for Solving Economic Dispatch Problems : A Postoptimization Analysis

This paper proposes a predictor-corrector primal-dual interior point method which introduces line search procedures IPLS in both the predictor and corrector steps. The Fibonacci search technique is used in the predictor step, while an Armijo line search is used in the corrector step. The method is developed for application to the economic dispatch ED problem studied in the field of power systems analysis. The theory of the method is examined for quadratic programming problems and involves the analysis of iterative schemes, computational implementation, and issues concerning the adaptation of the proposed algorithm to solve ED problems. Numerical results are presented, which demonstrate improvements and the efficiency of the IPLS method when compared to several other methods described in the literature. Finally, postoptimization analyses are performed for the solution of ED problems.


Introduction
Since its introduction in 1984, the projective transformation algorithm proposed by Karmarkar in 1 has proved to be a notable interior point method for solving linear programming problems LPPs .This pioneer study caused an upheaval in research activities in this area.Among all the variations of Karmarkar's original algorithm, the first to attract the attention of researchers was the one that uses a simple affine transformation to replace Karmarkar's original and highly complex projective transformation, enabling work with

The Economic Dispatch Problem
The economic dispatch problem ED is defined as an optimal allocation process of electricity demands among available generating units, where operational constraints must be satisfied while minimizing generation costs.Happ reported in 28 that, by 1920, several engineers had become aware of the economic allocation problem.According to the aforementioned author, one of the first methods employed to solve the ED problem was to request power from the most efficient unit the merit order loading method .This method was based on the following idea: the next incremental active power was to be supplied by the most efficient plant, until it reached its maximum operational point, and so on, successively.Although this method failed to minimize costs, it was employed until 1930, when the equal incremental cost criterion began to produce better results.
The idea behind the method of incremental costs is that the next incremental system load increase should be allocated to the unit with the lowest incremental cost, which is determined by measuring the derivative of the cost curve.Steinberg and Smith mathematically proved in 29 the equal incremental costs criterion, which was already being used empirically.Around 1931, this method had already become well established.In theory, the method described by Steinberg and Smith in 29 also ensures that if there are no active constraints at the point of optimal operation, the incremental costs of all units should be equal.This rule is still widely used by power system operators today.

Optimization Model for the Classical Economic Dispatch Problem
The economic dispatch ED problem is concerned with minimizing active power production costs while meeting the system demand and taking into account the operational limits of all generating units.The ED problem is mathematically described in where C T : total cost function for generating units, n: number of generating units, C i P i : cost of generation unit i without considering the valve-point effect , a i , b i , c i : coefficient of the cost function for generating unit i, P i : power output of generating unit i, D: total power demand, P min i : minimum power output limit for generating unit i, and P max i : maximum power output limit for generating unit i.
The objective function in 2.1 may also represent the so-called valve-point effects, as in 28 , which are associated with the opening of pressure valves at some specific operating points.In such cases, C i P i is mathematically described as in C i P i a i P 2 i b i P i c i e i sin f i P min i − P i .

2.2
The cost function C i P i described in 2.2 is continuous but not differentiable.Although it is more representative, function 2.2 makes ED a much more complex problem to solve due to its characteristics of nondifferentiability.
The literature describes several different methodologies to solve ED problems.In those studies, the methodology associated with evolutionary algorithms stands out, especially when issues related to nondifferentiability such as those described in 2.2 are involved.Evolutionary algorithms have been used to solve ED problems because they are able to find optimal solutions even when the objective function and/or constraints are not continuous or non-differentiable.Numerical problems related to evolutionary algorithms involve the inability to verify the optimality conditions associated with the solutions obtained and also the computational effort necessary to obtain the solutions, especially for large systems.
Optimal solutions to ED problems have also been investigated through traditional nonlinear programming methods, including interior point methods 6, 23, 30 .This paper proposes a predictor-corrector primal-dual interior point method for solving ED problems, which incorporates line search procedures in both the predictor and corrector steps.The results presented here demonstrate that the proposed method improves the solution of ED problems.The following section examines the application of the primal-dual interior point method to solve quadratic programming problems.

Solution Technique for the Quadratic Programming Problem
Quadratic programming problems QPPs represent a special class of nonlinear programming in which the objective function is quadratic and the constraints are linear 31 .The problem is expressed mathematically by

3.1
where and Q ∈ n×n is a diagonal matrix.Problem 3.1 is equivalent to 3.2 , where z ∈ n is a slack variable and r ∈ n is a surplus variable:

3.2
For μ > 0, it is possible to incorporate a logarithmic barrier function to the objective function of 3.2 and eliminate inequality constraints.This procedure results in the following nonlinear optimization problem: x − r l.

3.3
The Lagrangian function related to 3.3 is expressed by

3.4
where the dual variables associated with the three equality constraints in 3.3 are, respectively, w ∈ m , y ∈ n , s ∈ n .The optimal Karush-Kuhn-Tucker KKT conditions for problem 3.1 are depicted in where R, Z, S and Y are diagonal matrices whose diagonal elements are r i , z i , s i , and y i ; i 1, . . ., n, respectively; e 1, . . ., 1 T ; μ is a dual metric or adjustment parameter for the curve defined by the central trajectory path-following parameter .The set Ω 0 given in 3.6 is defined to simplify to notation.This set describes interior points for problem 3.2 and its corresponding dual problem: x, w, z, r, y, s ⊥ − Qx A T w s − y c, Ax b, x z u, x − r l, z, r, y, s > 0 .

3.6
Equations 3.4 , 3.5 , and 3.6 are considered in the following sections to perform the analysis of important issues concerning the proposed IPLS method, such as search directions, step sizes, stopping criteria, and updating of the barrier parameter.

Search Directions
In this section, the search directions used in the proposed IPLS method are investigated.In Section 3.2, the search directions for the predictor step are calculated while the search directions for the corrector step are determined in Section 3.3.The proposed strategy is a variant of the approach developed by Wu et al. in 23 .

Search Directions-Predictor Step
Let us suppose that, in an iteration k, a point h k satisfies the primal and dual feasibility conditions expressed by 3.5 .In this case, the definition of the new point h k 1 depends solely on the calculation of a search direction and on a step size in such a direction.Disregarding the step size in the iteration k 1, the new point h k 1 is defined by the following equation: Therefore, it is necessary to determine the direction of the movement d k to obtain the new point h k 1 .Following the steps in the Newton method applied to the nonlinear system 3.5 , d k can be obtained by solving for the system 3.8 , which is equivalent to 3.9 : where the residuals of 3.9 for the predictor step are expressed as in 14 15 From 3.13 and 3.14 , we obtain 3.17 and 3.18 , respectively: Isolating d k y and d k s in 3.15 and 3.16 , respectively, leads to 3.19 and 3.20 , as follows:

3.20
Combining the results found in 3.17 , 3.18 , 3.19 , and 3.20 with those given in 3.11 -3.16 , and considering 3.21 , yields the directions 3.22 and 3.23 : where Note that, after calculating 3.23 , the remaining components of the direction vector d k z , d k r , d k y , and d k s in 3.17 , 3.18 , 3.19 , and 3.20 , respectively, are easily calculated.Since matrix AθA T is symmetrical and positive definite in 3.22 considering that Q is symmetrical and positive definite , d k w can be determined by using the Cholesky decomposition.

Search Directions-Corrector Step
Analogously to the procedure carried out in Section 3.2, this section describes the calculation of the search direction for the corrector step d k , which is obtained by solving the linear system 3.25 , as follows: where F h k is obtained by considering second-order approximations in the residuals 3.10 calculated in the predictor step , so that 3.25 is equivalent to where The calculation of the residuals t k , g k , f k , o k , q k , v k , described in 3.10 , and of t k , g k , f k , o k , q k , v k , described in 3.27 , basically distinguishes the predictor and corrector steps in the proposed method.It is important to note that the corrector step procedure uses direction values d k x , d k z , d k y , and d k s , which have already been calculated in the predictor step, to redefine residuals q k and v k in 3.27 .Using 3.25 and 3.26 and following the same steps Mathematical Problems in Engineering taken to determine the directions of the predictor step, as seen in Section 3.2, the components of the direction vector d k can be calculated using the following: where and θ is defined in 3.21 .

Step Size
After calculating the search directions for the predictor and corrector steps, it is possible to move to a new point x k 1 , w k 1 , z k 1 , r k 1 , y k 1 , s k 1 , while ensuring that s k 1 > 0, z k 1 > 0, y k 1 > 0 and r k 1 > 0. In order to ensure nonnegativity constraints over the slack variables, the step to be taken in each direction in both the predictor and corrector steps must be controlled.The basics of this procedure are described in 31 and are discussed below.

Predictor Step Size
Considering the variables defined in the predictor step, the step size for primal and dual variables is calculated as described below:

3.40
The step size α P k for the primal variables is calculated through where, for 0 < α < 1, the step sizes α P , α Q , α LS are determined as follows.
i The step size for primal variables α P is obtained without violating the nonnegativity requirements of the primal variables: ii The step size α Q is the maximum step size possible without increasing the objective value:

3.43
iii The step size α LS is determined as in 3.44 from the Fibonacci line search strategy, which is briefly summarized as follows: where F μ is defined in 3.3 .
For notation simplicity, in the summary of the Fibonacci search, the following identities are defined: and d k z is defined in 3.17 .Only primal variables are considered in the Fibonacci search algorithm used by the IPLS method.Starting at a point x k x k , r k , z k T , the algorithm searches a new point x k 1 in direction d k x using the function F μ defined in 3.3 .The Fibonacci method calculates α LS so that the minimization of F μ is ensured.α LS is determined considering the Fibonacci sequence.The initial value of α LS is set taking into account the interval of uncertainty 0, 1 .
The step size α D k for the dual variables 3.38 -3.40 is calculated through 3.45 for 0 < α < 1: Mathematical Problems in Engineering

Corrector Step Size
In the corrector step, the calculation of step size is defined analogously, considering the variables from the following where while the step sizes β P , β Q , β LS are determined as follows.
i The step size for primal variables β P is obtained without violating the nonnegativity requirements of the primal variables: ii The step size β Q is the maximum step size possible without increasing the objective value:

3.54
iii The step size β LS in 3.55 is determined from the Armijo line search: are defined in 3.29 -3.31 , respectively.In the proposed IPLS method, the Armijo search is also performed in the direction defined by primal variables.Starting at a point x k x k , r k , z k , the method searches for a new point x k 1 , in direction d k x , so that the function F μ defined in 3.3 decreases.This search calculates β LS , thereby ensuring the reduction of F μ .To prevent oscillations in the iterative process, the initial choice for β k should not be too high, and to prevent the process from stopping prematurely, it should not be too low.Therefore, this value is generally adjusted to β 0 1.If 3.44 is not satisfied, β k is updated using the following sequence:

3.57
The following equation is used for the analysis of 3.56 : with μ k determined as in Section 3.6.
The step size for dual variables β D k is calculated by 3.59 , without violating the nonnegativity requirements of the dual variables:

3.59
The general principle used here to calculate α P k and β P k is to choose a step size that reduces the quadratic objective by a maximum amount without violating the nonnegativity requirements of the primal variables.
The basic idea for defining the line searches in both the predictor and corrector steps is to use a more accurate search for the predictor step which uses first order approximation to calculate the residuals and a simpler search, albeit more robust, for the corrector step which uses second-order approximation to calculate the residuals .Therefore, the Fibonacci search is used in the predictor step, since it is more accurate, and provides the minimum value for the objective function in the predefined direction; the Armijo search is used in the corrector step because it is simpler and more robust.

Stopping Rules
Interior point algorithms do not find exact solutions for linear or quadratic programming problems.Therefore, stopping rules are needed to decide when the solution obtained in a current iteration is sufficiently close to the optimal solution.In this study, the stopping rules are based on 32 .
Many algorithms consider that a good approximate solution is the one that presents sufficiently small values for primal and dual residuals t k , u k , and also for the dual metric μ k .Nevertheless, it is possible to use relative values for the metrics t k , u k , and μ k in order to Mathematical Problems in Engineering reduce scaling problems, as described in 32 .Typical stopping rules are shown in 3.60 -3.65 , as follows: i primal feasibility: ii dual feasibility: iii complementary slackness: where ε 1 , ε 2 , and ε 3 are sufficiently small positive numbers.Eventually, other criteria can be adopted according to the specific characteristics of each problem, as can be seen in 24, 32 .

Update of the Barrier Parameter
According to 32 , the barrier parameter is updated using an inner product that involves the primal variables r k and z k , and the dual variables s k and y k , respectively, as shown in the following equation 3.66 so that μ k is calculated using the following equation where the parameter σ is used to accelerate the convergence of the iterative process.This procedure for updating the barrier parameter proposed by Wright in 32 helps in the theoretical convergence proof and also in the complexity analysis of primal-dual methods.
Step 2 intermediate calculations-predictor .Calculate: t k ; g k ; f k ; o k ; q k ; v k using 3.10 , μ k using 3.67 and matrix θ using 3.21 .
Step 3 finding directions of translation-predictor .Calculate search directions and d k s for the predictor step using 3.17 -3.24 .
Step 4 computing step size-predictor .Calculate Fibonacci step sizes α P k and α D k using 3.41 -3.44 .
Step 5 moving to a new solution-predictor .Update x k 1 ; w k 1 ; s k 1 ; y k 1 ; z k 1 ;r k 1 obtained from the predictor step, according to 3.35 -3.44 .
Step 6 checking for optimality-predictor .If the criteria defined in 3.60 , 3.62 , and 3.64 are satisfied, stop.The solution is optimal.Otherwise, go on to the following step.
Step 7 intermediate calculations-corrector .Calculate: t k ; g k ; f k ; o k ; q k ; v k using 3.27 , μ k using 3.67 and matrix θ using 3.21 .
Step 8 finding directions of translation-corrector .Calculate and d k s , for the corrector step using 3.28 -3.33 .
Step 9 computing step size-corrector .Calculate the Armijo step size β P k and β D k using 3.52 -3.58 .
Step 11 checking for optimality-corrector .If the criteria defined in 3.61 , 3.63 , and 3.65 are satisfied, stop.The solution is optimal.Otherwise, go on to the next step.
Step 12. Adjust k k 1 and return to Step 2.

The predictor Steps 2 through
Step 6 are held in odd iterations, while the corrector Steps 7 through Step 11 are performed in even iterations.The next section describes numerical simulations involving the application of the proposed method to ED problems.

Application of the Proposed Algorithm to Solve ED Problems
In this section, the proposed method is applied to solve three power systems with 3, 6, and 13 generators, respectively.Tables 1, 3, and 5 show the data related to the generating units of the power systems.These data were extracted from 26, 27 .The tables list the coefficients of the cost function for each generating unit, as described in 2.1 , and also minimum and maximum power output limits.
Tables 2 and 3 present the results of the application of the proposed method to solve the power system with 3 generators, while Tables 5 and 6 list the results for the system with 6 generators, and Tables 8 and 9 depict the results for the system with 13 generators.These tables compare the solutions obtained by the proposed IPLS method against those calculated by the following methods: the predictor-corrector primal-dual PCPD method described in 6 , the hybrid genetic algorithm HGA , the coevolutionary genetic algorithm COEGA , the hybrid atavistic genetic algorithm HAGA and the cultural algorithm CA .The solutions determined by the HGA and COEGA methods are available in 27 , while the solutions obtained with the CA method are described in 26 , and those obtained with HAGA method are given in 25 , but only for the system comprising 13 generators.
The comparison between the IPLS method and the evolutionary approaches cited above were introduced in this section because these are traditional methods used in power system to solve ED problems, especially when a more general nondifferentiable objective function as shown in 2.2 is used.However, it is important to highlight that these evolutionary approaches are heuristic procedures, which provide only approximate solutions to ED problems.When ED is formulated as a quadratic problem, it can be solved by means of exact methods, such as the IPLS and PCPD method.
Therefore, to better evaluate the performance of proposed IPLS method, this method has been compared to the PCPD method in Table 10.The results in Table 10 have the main purpose to show the reduction in the computational effort when the IPLS method is compared to the PCPD method.Both methods were implemented using Borland Pascal 7.0 programming language.

Power System with 3 Generators
The main characteristics of the system containing 3 generators are described in Table 1.Parameters a i , b i , c i stand for the coefficients of the cost functions for the generators, while P min i , P max i represent minimum and maximum power output capabilities, respectively, for generating unit i.

4.1
The parameters related to the system's total demand are set at D 850 MW and the active power losses are neglected.The values adopted for ε 1 , ε 2 , and ε 3 are 10 −8 .
Table 2 shows the active power output calculated by the methods.The results for the HAGA algorithm are not presented in this case study.Table 3 compares optimal values for the objective functions obtained by each method.From the results presented in this table, it is clear that the dispatches calculated by all the types of genetic algorithms cannot reach the global optimum dispatch attained by the interior point methods PCPD and IPLS.As Table 3 indicates, the cost calculated by the IPLS and PCPD methods is lower than that obtained by the HGA and COEGA methods.As already discussed, this is an expected result, since the evolutionary approaches provide only approximate solution to the problem.

Power System with 6 Generators
Table 4 describes the main characteristics of the system containing 6 generators.
As in the previous case study, the parameters a i , b i , c i stand for the coefficients of the cost functions for the generating unit I, while P min i , P max i represent minimum and maximum power output capabilities, respectively, for the generating unit i.
Table 5 shows the active power dispatch calculated by the methods.The solutions for COEGA and HAGA algorithms were not presented by their authors.Table 6 compares optimal values for the objective function obtained by each method.Again, as expected, the costs calculated by the IPLS and PCPD methods are lower than those obtained by the HGA method.

4.3
The parameters related to the system's total demand are set at D 2520 MW and the active power losses are neglected.The values adopted for ε 1 , ε 2 , and ε 3 are 10 −8 .both methods tested.Therefore, the efficiency of the methods is compared only in terms of number of iterations, as described in Table 10.
The basic difference between the PCPD and IPLS algorithms is that line searches are incorporated in the corrector and predictor steps of the IPLS method.Therefore, these results demonstrate that the introduction of line searches in the IPLS method greatly reduces the number of iterations required to solve ED problems.Since line searches represent only a minor additional computational effort in the overall process of finding a solution, the computational effort for solving ED problems is significantly improved by the proposed IPLS method.

Postoptimization Analysis: Incremental Costs
This section conducts a postoptimization analysis to evaluate the optimality conditions of the solutions obtained in the previous section for the systems in question.The purpose of this analysis is to evaluate the complete picture concerning optimal solutions for ED problems.The classical ED problem 2.1 can be rewritten as follows: where Problem 4.4 is subsequently used for the postoptimization analysis of the solutions obtained by the systems studied in the previous section.

The Lagrangian Function and KKT Conditions
Considering ED as redefined in 4.4 , one finds the following associated Lagrangian function:

4.10
In the literature on power systems, the Lagrange multiplier w is called energy price i.e., the price of 1 MWh , while λ i P i is called incremental or marginal cost of generating unit i.Using expression 4.9 , the energy price can be calculated by w λ i P i − s i y i ; i 1, . . ., n. 4.11

Energy Prices and Incremental Costs Using KKT Conditions
If the optimal solution provided by the IPLS method satisfies the KKT conditions, then there are four possible combinations for analysing prices and incremental costs, which are examined below.

No Active Constraint
In this case there is no active constraint in the optimal solution of problem 4.4 , so that s * i 0, y * i 0; i 1, . . ., n.Thus, from 4.11 , one reaches w λ * i P * i ; i 1, . . ., n, 4.12 that is, the energy price is equal to the incremental marginal costs for all the generating units.This situation corresponds to the rule commonly used by power system operators, which states that all marginal costs are equal.It is important to emphasize that this rule is valid only for this case and should, therefore, be used cautiously.

Active Constraints for Maximum Power Output Limit
In this case, some generating units are dispatched in their maximum power output capabilities, so that the constraints associated with maximum power output capabilities are active in the optimal solution of 4.4 .To analyse this case, let the set Ω max be defined with the indices of the generating units that have been optimally dispatched in their maximum power output.In this situation, s * i 0 for i 1, . . ., n, y * j > 0 for j ∈ Ω max , and y * j 0 for j / ∈ Ω max .Starting from 4.11 , one reaches w λ * j P * j y * j ; j ∈ Ω max , w λ * j P * j ; j / ∈ Ω max .

4.13
Based on 4.13 , it is possible to show that the marginal costs for the group of generating units that belong to Ω max are lower than the marginal costs for the group that does not belong to Ω max , since y * j > 0; j ∈ Ω max .

Active Constraints for Minimum Power Output Limit
In this case, some generating units are dispatched in their minimum power output capabilities, so that the constraints associated with minimum power output capabilities are active in the optimal solution of 4.4 .To analyse this case, let the set Ω min be defined with the indices of the generating units that have been optimally dispatched in their minimum power output.In this situation, y * i 0 for i 1, . . ., n, s * j > 0 for j ∈ Ω min , and s * j 0 for j / ∈ Ω min .Thus, using 4.11 , one reaches w λ * j P * j − s * j ; j ∈ Ω min , w λ * j P * j ; j / ∈ Ω min .

4.14
Based on 4.14 , it can be shown that the marginal costs for the group of generating units that belong to Ω min are higher than the marginal costs for the group that does not belong to Ω min , since s * j > 0; j ∈ Ω min .

Active Constraints for Both Minimum and Maximum Power Output Capabilities
In this case, some generating units are dispatched in their upper limits, while others are dispatched in their lower limits, or at some other operational point between their upper and lower limits.Obviously, the same unit cannot simultaneously assume lower and upper limits.The comments made in the preceding sections concerning marginal costs for the generating units that are dispatched in their lower or upper limits also apply to this case.lower or upper limit, so there is no active constraint in the optimal solution.This situation coincides with the one described in Section 4.7.1.The energy price, marginal costs, and Lagrange multipliers are shown in Table 11 for this case.As described in Section 4.7.1, the values obtained by the IPLS method for s * i and y * i ; i 1, . . ., n are all zero.Also, the energy price is equivalent to the marginal costs, which are all equal, as expected.These results are in accordance with the theory described in Section 4.7.1.

Case 2: Power System with 6 Generating Units
Table 5 shows the optimal dispatch, P * i , i 1, . . ., n, calculated by the IPLS method for the system with 6 generating units.As can be seen in this table, the generating unit 2 has reached its lower limit, so the optimal solution has one active constraint.This situation coincides with the one described in Section 4.7.3.Table 12 lists the energy price, marginal costs, and Lagrange multipliers calculated by the IPLS method for this case.As described in Section 4.7.3, the values obtained by the IPLS method for s * i and y * i ; i 1, . . ., n are all zero except for the value of s * 2 , which is associated with the constraint for lower output capabilities at unit 2. Except for unit 2, all the marginal costs are equal to the energy price.These results are in accordance with the theory described Section 4.7.3 and also with 4.14 .

Case 3: Power System with 13 Generating Units
Table 8 shows the optimal dispatch, P * i , i 1, . . ., n, calculated by the IPLS method for the system with 13 generating units.Table 13 lists the energy price, marginal costs, and Lagrange multipliers calculated by the IPLS method.As the latter table indicates, generated units 1, 2, and 3 have reached their upper limits, and generated units 10, 11, 12, and 13 have reached their lower limits, indicating that there are 6 active constraints in the optimal solution calculated by the IPLS method.This situation coincides with the one described in Section 4.7.The analysed results indicate that the dispatches P * i , i 1, . . ., n calculated, respectively, in Tables 2, 5, and 8 and the results for λ * i , w, s * i , and y * i determined, respectively, in Tables 11,12, and 13, satisfy the complementary slackness conditions set forth in 4.7 and 4.8 .

Conclusions
This paper proposes a predictor-corrector primal-dual interior point method which introduces line search procedures IPLS in both the predictor and corrector steps.The Fibonacci search is used in the predictor step, while an Armijo search is used in the corrector step.The method is developed for application to the economic dispatch ED problem, which is an example of a quadratic programming problem studied in the field of power systems analysis.ED problems have already been solved through primal-dual interior point methods, although the strategy adopted here to solve the problem has not yet been tested.
The proposed algorithm is applied to solve ED problems in case studies involving systems with 3, 6, and 13 generating units.The results provided by the proposed IPLS method are compared to those provided by several other methods described in the literature, such as the predictor-corrector primal-dual PCPD interior point method, hybrid genetic algorithm, coevolutionary genetic algorithm, hybrid atavistic genetic algorithm, and cultural algorithm.The dispatches calculated by all the types of genetic algorithms could not reach the global optimal dispatch attained by the interior point methods PCPD and IPLS.
The computational effort of the interior point methods tested PCPD and IPLS was evaluated and measured in terms of the number of iterations required to find the optimal solution.The results indicate that the introduction of line searches in the IPLS method considerably reduces the number of iterations required for solving ED problems.Line searches pose represent only a minor additional computational effort in the overall process of finding a solution; hence, the computational effort for solving ED problems is also greatly improved by the proposed IPLS method.
A theory for performing a postoptimization analysis was also analysed.This theory highlights the relation among energy prices, incremental generation costs, and other Lagrange multipliers.This mathematical relation is used to validate optimal solutions for ED problems calculated by the IPLS method.
Further research could involve the representation of the so-called valve point loading in the objective function of ED problems.In that case, the nondifferentiability of the objective function should be treated appropriately.

Table 1 :
Characteristics of the system with 3 generators.

Table 2 :
Comparison of the power generation outputs obtained for the system with 3 generators.

Table 3 :
Characteristics of the system with 3 generators.

Table 4 :
Characteristics of the system with 6 generators.

Table 6 :
Values of the objective function for the system with 6 generators.

Table 7 :
Characteristics of the system with 13 generators.

Table 8 :
Comparison of the power generation outputs obtained for the system with 13 generators.

Table 10 :
Number of iterations for the PCPD and IPLS methods.
, s i ∈ , and y i ∈ , i 1, . . ., n, are the associated Lagrange multipliers.The KKT conditions associated with the Lagrangian function 4.5 are expressed through i ; i 1, . . ., n.

Table 2
describes the optimal dispatch, P * i , i 1, . . ., n, calculated by the IPLS method for the system with 3 generating units.As this table indicates, no generating unit has reached its

Table 11 :
Incremental analysis for the 3-generator system.

Table 12 :
Incremental analysis of the 6-generator system.
4. As described previously, the values obtained by the IPLS method for s * i and y * i ; i 1, . . ., n are zero except for the values of y * 1 , y * 2 , y * 3 , which are associated with upper output capabilities, and also for s * 10 , s * 11 , s * 12 , s * 13 , which are associated with lower output

Table 13 :
Incremental analysis for the 13-generator system.The marginal costs for all the remaining units 4, 5, 6, 7, 8, and 9 are equal to the energy price.These results are in accordance with the theory described in Sections 4.7.2 and 4.7.3 and also with 4.13 and 4.14 .As expected, in all the cases analysed in Tables11, 12, and 13, the value of the generation cost w λ * i − s * i y * i was determined univocally, which proves the equal incremental cost criterion 29 .