`Mathematical Problems in EngineeringVolume 2012, Article ID 376546, 26 pageshttp://dx.doi.org/10.1155/2012/376546`
Research Article

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

Received 9 November 2011; Accepted 3 April 2012

Copyright © 2012 Antonio Roberto Balbo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

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.

#### 1. 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 LPP in its standard form. The affine algorithm was first introduced in [2] by Dikin, a Soviet mathematician. Later, in 1986, the work was independently rediscovered by Barnes in [3] and by Vanderbei et al. in [4]. They proposed using the primal-affine algorithm to solve LPP in standard form and also presented proof of the algorithm’s convergence. A similar algorithm, called dual-affine, was developed and implemented by Adler et al. [5] to solve the LPP in the form of inequality. Compared to the relatively cumbersome projective transformation, the implementation of the primal-affine and dual-affine algorithms was simpler because of its direct relationship with the LPP. These two algorithms produced promising results when applied to large problems [6], although the theoretical proof of polynomial time complexity was not obtained from the affine transformation. Megiddo and Shub’s study in [7] indicated that the trajectory that leads to the optimal solution provided by affine algorithms depends on the initial solution. A poor initial solution, which is close to a viable domain vertex, could result in an investigation that covers all vertex problems.

Nevertheless, the polynomial time complexity of primal-affine and dual-affine algorithms can be re-established by incorporating a logarithmic barrier function into the objective function of the original LPP. The purpose of this procedure is to solve the problem pointed out by Megiddo, that is, to prevent an interior solution from becoming “trapped” at the border of the problem (possibly a vertex). This procedure also provides proof of the complexity of the method. The idea of using the logarithmic barrier function method for convex programming problems was developed by Fiacco and McCormick in [8, 9], based on the method proposed by Frisch in [10]. After the introduction of Karmarkar’s algorithm in 1984, the logarithmic barrier function method was reconsidered to solve linear programming problems. Gill et al. used this method in [11] to develop a projected Newton barrier method and demonstrated its equivalence with Kamarkar’s projective algorithm. The methods proposed, among others, by Ye in [12], Renegar in [13], Vaidya in [14], and Megiddo in [15], as well as those of a central trajectory—called path-following methods, which were proposed by Gonzaga in [16, 17] and Monteiro and Adler in [18], use the objective function augmented by the logarithmic barrier function. A third variation, the so-called primal-dual interior point algorithm, was introduced by Monteiro et al. in [19] and by Kojima et al. in [20]. This algorithm explores a potential primal-dual function, a variation of the logarithmic barrier function, called the potential function. The polynomial time complexity theory was successfully demonstrated by Kojima et al. in [20] and Monteiro et al. in [19] based on Megiddo’s work, which provided a theoretical analysis for the logarithmic barrier method and proposed the primal-dual approach.

The predictor-corrector procedure initially defined by Mehrotra and Sun in [21] and implemented in by Lustig et al. in [22] explored variant directions of the primal-affine interior point methods in the predictor step. In the corrector step, the point previously obtained (in the predictor step) was “centralized” to exploit the potential function related to the logarithmic barrier function. This procedure significantly improved the performance of the primal-dual interior point methods.

This strategy was reviewed and modified by Wu et al. in [23] and successfully applied to solve optimal power flow problems. Wu used the logarithmic barrier penalization, the Newton method, and first-order approximations in the predictor step to determine search directions and the approximate solution. Second-order approximations were considered in the corrector step to refine solutions obtained in the predictor step.

The methods related to the primal-dual interior point methodology, especially those proposed by Kojima et al. in [20], Monteiro and Adler in [18], and Monteiro et al. in [19], which were broadly investigated by Fang and Puthenpura in [24], have been explored in this past decade to solve linear, nonlinear, and integer mathematical programming problems in several fields of research.

This paper proposes a primal-dual interior point method that uses the predictor-corrector strategy described by Lustig et al. in [22] and Wu et al. in [23] and incorporates line search procedures in both the predictor and corrector steps. Theoretical aspects as well as iterative schemes and computational implementation are investigated. A Fibonacci line search procedure is carried out in the predictor step and an Armijo line search is used in the corrector step to calculate the step sizes while taking into account constraints in the variables of the problem. The line search procedures adopted are aimed at improving the overall convergence of the proposed method for solving quadratic programming problems. The method is applied to solve the economic dispatch (ED) problem, a classical quadratic problem studied in the field of electrical power systems. The results obtained with the proposed method are compared with those obtained by several others described in the literature for solving ED problems, such as the primal-dual method described in [6], evolutionary algorithms found in [2527], genetic and coevolutionary genetic algorithms described by Samed in [27], the cultural algorithm described by Rodrigues in [26], and a hybrid atavistic genetic algorithm given in [25]. This comparative investigation demonstrates the efficiency of the proposed IPLS method.

This paper is organized as follows: Section 2 presents the ED problem; Section 3 develops the theory of the proposed IPLS method and presents its algorithm; Section 4 describes a computational implementation of the method, applying it to solve ED problems. This section also includes a postoptimization analysis. Finally, conclusions are drawn in Section 5.

#### 2. 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.

##### 2.1. 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 : total cost function for generating units, : number of generating units, : cost of generation unit (without considering the valve-point effect), : coefficient of the cost function for generating unit , : power output of generating unit , : total power demand, : minimum power output limit for generating unit , and : maximum power output limit for generating unit .

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, is mathematically described as in

The cost function 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.

#### 3. 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 where so that ,  , and is a diagonal matrix. Problem (3.1) is equivalent to (3.2), where is a slack variable and is a surplus variable:

For , 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:

The Lagrangian function related to (3.3) is expressed by where the dual variables associated with the three equality constraints in (3.3) are, respectively, . The optimal Karush-Kuhn-Tucker (KKT) conditions for problem (3.1) are depicted in where and are diagonal matrices whose diagonal elements are , and ;, respectively; ; is a dual metric or adjustment parameter for the curve defined by the central trajectory (path-following parameter). The set given in (3.6) is defined to simplify to notation. This set describes interior points for problem (3.2) and its corresponding dual problem:

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.

##### 3.1. 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].

##### 3.2. Search Directions—Predictor Step

Let us suppose that, in an iteration , a point satisfies the primal and dual feasibility conditions expressed by (3.5). In this case, the definition of the new point 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 , the new point is defined by the following equation:

Therefore, it is necessary to determine the direction of the movement to obtain the new point . Following the steps in the Newton method applied to the nonlinear system (3.5), 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

The directions to be determined in the predictor step use the residuals defined in (3.10), resulting in (3.11)–(3.16), as follows:

From (3.13) and (3.14), we obtain (3.17) and (3.18), respectively:

Isolating and in (3.15) and (3.16), respectively, leads to (3.19) and (3.20), as follows:

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 , and in (3.17), (3.18), (3.19), and (3.20), respectively, are easily calculated. Since matrix is symmetrical and positive definite in (3.22) (considering that is symmetrical and positive definite), can be determined by using the Cholesky decomposition.

##### 3.3. 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, which is obtained by solving the linear system (3.25), as follows: where 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, , and are diagonal matrices whose diagonal components are , , and , , respectively.

The calculation of the residuals , described in (3.10), and of , 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 , , , and , which have already been calculated in the predictor step, to redefine residuals and in (3.27). Using (3.25) and (3.26) and following the same steps taken to determine the directions of the predictor step, as seen in Section 3.2, the components of the direction vector can be calculated using the following: where and is defined in (3.21).

##### 3.4. Step Size

After calculating the search directions for the predictor and corrector steps, it is possible to move to a new point , while ensuring that , , and . 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.

###### 3.4.1. Predictor Step Size

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

The step size for the primal variables is calculated through where, for , the step sizes are determined as follows.(i)The step size for primal variables is obtained without violating the nonnegativity requirements of the primal variables:(ii)The step size is the maximum step size possible without increasing the objective value:(iii)The step size is determined as in (3.44) from the Fibonacci line search strategy, which is briefly summarized as follows: where is defined in (3.3).

For notation simplicity, in the summary of the Fibonacci search, the following identities are defined: , and ; so that is defined in (3.23), is defined in (3.18), and is defined in (3.17). Only primal variables are considered in the Fibonacci search algorithm used by the IPLS method. Starting at a point , the algorithm searches a new point in direction using the function defined in (3.3). The Fibonacci method calculates so that the minimization of is ensured. is determined considering the Fibonacci sequence. The initial value of is set taking into account the interval of uncertainty .

The step size for the dual variables (3.38)–(3.40) is calculated through (3.45) for :

###### 3.4.2. 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 are determined as follows.(i)The step size for primal variables is obtained without violating the nonnegativity requirements of the primal variables:(ii)The step size is the maximum step size possible without increasing the objective value:(iii)The step size in (3.55) is determined from the Armijo line search: so that where ; ; so that 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 , the method searches for a new point , in direction , so that the function defined in (3.3) decreases. This search calculates , thereby ensuring the reduction of . To prevent oscillations in the iterative process, the initial choice for 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 . If (3.44) is not satisfied, is updated using the following sequence:

The following equation is used for the analysis of (3.56): with determined as in Section 3.6.

The step size for dual variables is calculated by (3.59), without violating the nonnegativity requirements of the dual variables:

The general principle used here to calculate and 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.

##### 3.5. 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 , and also for the dual metric . Nevertheless, it is possible to use relative values for the metrics , and in order to 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 , , and 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].

##### 3.6. Update of the Barrier Parameter

According to [32], the barrier parameter is updated using an inner product that involves the primal variables and , and the dual variables and , respectively, as shown in the following equation so that 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.

##### 3.7. Algorithm for the Proposed IPLS Method

Step 1 (initialisation). Adjust . Choose an arbitrary point: and choose , , and as sufficiently small positive numbers.

Step 2 (intermediate calculations—predictor). Calculate: using (3.10), using (3.67) and matrix using (3.21).

Step 3 (finding directions of translation—predictor). Calculate search directions , , , , , and for the predictor step using (3.17)–(3.24).

Step 4 (computing step size—predictor). Calculate Fibonacci step sizes and using (3.41)–(3.44).

Step 5 (moving to a new solution—predictor). Update ; ; ; ; ; 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: ; ; ; ; ; using (3.27), using (3.67) and matrix using (3.21).

Step 8 (finding directions of translation—corrector). Calculate , , , , , and , for the corrector step using (3.28)–(3.33).

Step 9 (computing step size—corrector). Calculate the Armijo step size and using (3.52)–(3.58).

Step 10 (moving to a new solution—corrector). Update ; ; ; ; ; using (3.46)–(3.48).

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.

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.

#### 4. 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.

Table 1: Characteristics of the system with 3 generators.

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.

Table 2: Comparison of the power generation outputs obtained for the system with 3 generators.
Table 3: Characteristics of the system with 3 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.

##### 4.1. Power System with 3 Generators

The main characteristics of the system containing 3 generators are described in Table 1. Parameters stand for the coefficients of the cost functions for the generators, while represent minimum and maximum power output capabilities, respectively, for generating unit .

The following values are adopted to initialise the method:

The parameters related to the system’s total demand are set at and the active power losses are neglected. The values adopted for ,, and are .

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.

##### 4.2. Power System with 6 Generators

Table 4 describes the main characteristics of the system containing 6 generators.

Table 4: Characteristics of the system with 6 generators.
Table 5: Comparison of the power generation outputs obtained for the system with 6 generators.
Table 6: Values of the objective function for the system with 6 generators.

As in the previous case study, the parameters stand for the coefficients of the cost functions for the generating unit , while represent minimum and maximum power output capabilities, respectively, for the generating unit .

The following values are adopted to initialise the method:

The parameters related to the system’s total demand are set at and the active power losses are neglected. The values adopted for ,, and are .

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. Power System with 13 Generators

Table 7 describes the main characteristics of the system containing 13 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 9: Values of the objective function for the system with 13 generators.
Table 10: Number of iterations for the PCPD and IPLS methods.

The parameters in this table are analogous to those described in the preceding case studies. The following values are adopted to initialise the method:

The parameters related to the system’s total demand are set at and the active power losses are neglected. The values adopted for ,, and are .

Table 8 lists the power generation outputs obtained by the methods. The values of the dispatch calculated by the HAGA are not given by [25], who provided only the optimal value for the objective function. Once more, as Table 9 indicates, the total cost calculated by the PCPD and IPLS methods is lower than that calculated by all the others, although CA and HAGA approaches get close to the optimal solution point.

##### 4.4. Performance of the Interior Point Methods Tested

This section evaluates the computational effort of the interior point methods tested here (PCPD and IPLS), measuring it in terms of the number of iterations required to obtain the optimal solution. To this end, Table 10 shows the number of iterations obtained by the solution algorithms of each method for the previously studied systems. It is important to highlight that the same parameter settings were used for the two methods.

Computational tests were carried out in an Intel Corel Quad Q9550, with 3.5 GB of RAM memory, in order to calculate and compare the CPU times between the proposed IPLS method and the PCPD method. For such a purpose, we utilized the specific unit GET TIME from Borland Pascal 7.0. The precision of this unit is up to milliseconds. The computational times calculated by this unit for the power systems tested (which include the systems with 3, 6, and 13 generators) were all null for both, the proposed method, and the PCPD method. This occurred due to the efficiency of the machine processor and also due to the efficiency of 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.

##### 4.5. 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 , and = = .

Problem (4.4) is subsequently used for the postoptimization analysis of the solutions obtained by the systems studied in the previous section.

##### 4.6. The Lagrangian Function and KKT Conditions

Considering ED as redefined in (4.4), one finds the following associated Lagrangian function: where ,  , and ,  , are the associated Lagrange multipliers. The KKT conditions associated with the Lagrangian function (4.5) are expressed through where

In the literature on power systems, the Lagrange multiplier is called energy price (i.e., the price of 1 MWh), while is called incremental or marginal cost of generating unit . Using expression (4.9), the energy price can be calculated by

##### 4.7. 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.

###### 4.7.1. No Active Constraint

In this case there is no active constraint in the optimal solution of problem (4.4), so that, ; . Thus, from (4.11), one reaches 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.

###### 4.7.2. 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 be defined with the indices of the generating units that have been optimally dispatched in their maximum power output. In this situation, for , for , and for . Starting from (4.11), one reaches

Based on (4.13), it is possible to show that the marginal costs for the group of generating units that belong to are lower than the marginal costs for the group that does not belong to , since .

###### 4.7.3. 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 be defined with the indices of the generating units that have been optimally dispatched in their minimum power output. In this situation, for , for , and for . Thus, using (4.11), one reaches

Based on (4.14), it can be shown that the marginal costs for the group of generating units that belong to are higher than the marginal costs for the group that does not belong to , since .

###### 4.7.4. 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.

##### 4.8. Results of Incremental Analyses for the Systems under Study
###### 4.8.1. Case 1: Power System with 3 Generating Units

Table 2 describes the optimal dispatch, , calculated by the IPLS method for the system with 3 generating units. As this table indicates, no generating unit has reached its 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 and ; 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.

Table 11: Incremental analysis for the 3-generator system.
###### 4.8.2. Case 2: Power System with 6 Generating Units

Table 5 shows the optimal dispatch, , 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 and ; are all zero except for the value of , 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).

Table 12: Incremental analysis of the 6-generator system.
###### 4.8.3. Case 3: Power System with 13 Generating Units

Table 8 shows the optimal dispatch, , 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.4. As described previously, the values obtained by the IPLS method for and; are zero except for the values of , which are associated with upper output capabilities, and also for , which are associated with lower output capabilities. 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).

Table 13: Incremental analysis for the 13-generator system.

As expected, in all the cases analysed in Tables 11, 12, and 13, the value of the generation cost was determined univocally, which proves the equal incremental cost criterion [29].

The analysed results indicate that the dispatches calculated, respectively, in Tables 2, 5, and 8 and the results for , , , and determined, respectively, in Tables 11, 12, and 13, satisfy the complementary slackness conditions set forth in (4.7) and (4.8).

#### 5. 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.

#### Acknowledgment

The authors are grateful for the financial support of FAPESP and CNPq (Brazil).

#### References

1. N. Karmarkar, “A new polynomial-time algorithm for linear programming,” Combinatorica, vol. 4, no. 4, pp. 373–395, 1984.
2. I. I. Dikin, “Iterative solution of problems of linear and quadratic programming,” Doklady Akademii Nauk SSSR, vol. 174, pp. 747–748, 1967 (Russian).
3. E. R. Barnes, “A variation on Karmarkar's algorithm for solving linear programming problems,” Mathematical Programming, vol. 36, no. 2, pp. 174–182, 1986.
4. R. J. Vanderbei, M. S. Meketon, and B. A. Freedman, “A modification of Karmarkar's linear programming algorithm,” Algorithmica, vol. 1, no. 4, pp. 395–407, 1986.
5. I. Adler, M. G. C. Resende, G. Veiga, and N. Karmarkar, “An implementation of Karmarkar's algorithm for linear programming,” Mathematical Programming, vol. 44, no. 1–3, pp. 297–335, 1989.
6. A. R. Balbo, M. A. S. Souza, and E. C. Baptista, “Métodos primal-dual de pontos interiores aplicados à resolução de problemas de despacho econômico: sobre a influência da solução inicial,” in Proceedings of the Simpósio Brasileiro de Pesquisa Operacional, João Pessoa—Pb. Anais do XL SBPO (XL SBPO '08), pp. 2074–2085, ILTC, Rio de Janeiro, RJ, Brazil, 2008.
7. N. Megiddo and M. Shub, “Boundary behavior of interior point algorithms in linear programming,” Mathematics of Operations Research, vol. 14, no. 1, pp. 97–146, 1989.
8. A. V. Fiacco and G. P. McCormick, Nonlinear Programming: Sequential Unconstrained Minimization Techniques, John Wiley & Sons, New York, NY, USA, 1968.
9. A. V. Fiacco and G. P. McCormick, Nonlinear Programming, vol. 4 of Classics in Applied Mathematics, SIAM, Philadelphia, Pa, USA, 2nd edition, 1990.
10. K. R. Frisch, The Logarithmic Potential Method of Convex Programming, University Institute of Economics, Oslo, Norway, 1955.
11. P. E. Gill, W. Murray, M. A. Saunders, J. A. Tomlin, and M. H. Wright, “On projected Newton barrier methods for linear programming and an equivalence to Karmarkar's projective method,” Mathematical Programming, vol. 36, no. 2, pp. 183–209, 1986.
12. Y. Ye, “An $O\left({n}^{3}L\right)$ potential reduction algorithm for linear programming,” in Contemporary Mathematics, vol. 114, pp. 91–107, American Mathematical Society, Providence, RI, USA, 1990.
13. J. Renegar, “A polynomial-time algorithm, based on Newton's method, for linear programming,” Mathematical Programming, vol. 40, no. 1, pp. 59–93, 1988.
14. P. M. Vaidya, “An algorithm for linear programming which requires $O\left(\left(\left(m+n\right){n}^{2}+{\left(m+n\right)}^{1.5}n\right)L\right)$ arithmetic operations,” Mathematical Programming, vol. 47, no. 2, pp. 175–201, 1990.
15. N. Megiddo, “On the complexity of linear programming,” in Advances in Economic Theory, T. Bewley, Ed., pp. 225–268, Cambridge University Press, Cambridge, UK, 1989.
16. C. C. Gonzaga, “An algorithm for solving linear programming problems in $O\left({n}^{3}L\right)$ operations,” in Progress in Mathematical Programming: Interior-Point and Related Methods, N. Megiddo, Ed., pp. 1–28, Springer, New York, NY, USA, 1989.
17. C. C. Gonzaga, “Polynomial affine algorithms for linear programming,” Mathematical Programming, vol. 49, no. 1, pp. 7–21, 1990.
18. R. D. C. Monteiro and I. Adler, “Interior path following primal-dual algorithms. I. Linear programming,” Mathematical Programming, vol. 44, no. 1, pp. 27–41, 1989.
19. R. D. C. Monteiro, I. Adler, and M. G. C. Resende, “A polynomial-time primal-dual affine scaling algorithm for linear and convex quadratic programming and its power series extension,” Mathematics of Operations Research, vol. 15, no. 2, pp. 191–214, 1990.
20. M. Kojima, S. Mizuno, and A. Yoshise, “A primal-dual interior point algorithm for linear programming,” in Progress in mathematical programming: Interior-Point and Related Methods, N. Megiddo, Ed., pp. 29–47, Springer, New York, NY, USA, 1989.
21. S. Mehrotra and J. Sun, “An algorithm for convex quadratic programming that requires $O\left({n}^{3.5}L\right)$ arithmetic operations,” Mathematics of Operations Research, vol. 15, no. 2, pp. 342–363, 1990.
22. I. J. Lustig, R. E. Marsten, and D. F. Shanno, “On implementing Mehrotra's predictor-corrector interior-point method for linear programming,” SIAM Journal on Optimization, vol. 2, no. 3, pp. 435–449, 1992.
23. Y. C. Wu, A. S. Debs, and R. E. Marsten, “Direct nonlinear predictor-corrector primal-dual interior point algorithm for optimal power flows,” IEEE Transactions on Power Systems, vol. 9, no. 2, pp. 876–883, 1994.
24. S. C. Fang and S. Puthenpura, Linear Optimization and Extensions: Theory and Algorithms, Prentice-Hall, Englewood Cliffs, NJ, USA, 1993.
25. J. O. Kim, D. J. Shin, J. N. Park, and C. Singh, “Atavistic genetic algorithm for economic dispatch with valve point effect,” Electric Power Systems Research, vol. 62, no. 3, pp. 201–207, 2002.
26. N. M. Rodrigues, Um algoritmo cultural para problemas de despacho de energia elétrica, Dissertação de Mestrado, Universidade Estadual de Maringá, 2007.
27. M. M. A. Samed, Um algoritmo genético Hibrido co-evolutivo para resolver problemas de despacho, Tese de Doutorado, UEM, Departamento de Engenharia Química, 2004.
28. H. H. Happ, “Optimal power dispatch—a comprehensive survey,” IEEE Transactions on Power Apparatus and Systems, vol. 96, no. 3, pp. 841–854, 1977.
29. M. J. C. Steinberg and T. H. Smith, Economic Loading of Power Plants and Electric Systems, MacGraw-Hill, 1943.
30. A. R. L. Oliveira, S. Soares, and L. Nepomuceno, “Optimal active power dispatch combining network flow and interior point approaches,” IEEE Transactions on Power Systems, vol. 18, no. 4, pp. 1235–1240, 2003.
31. M. S. Bazaraa, H. D. Sherali, and C. M. Shetty, Non Linear Programming-Theory and Algorithm, John Wiley & Sons, New York, NY, USA, 2nd edition, 1993.
32. S. J. Wright, Primal-Dual Interior-Point Methods, SIAM, Philadelphia, Pa, USA, 1997.