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

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 [25–27], 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 Minimize𝐶_𝑇=^𝑛_𝑖=1𝐶_𝑖𝑃_𝑖=𝑎_𝑖𝑃²_𝑖+𝑏_𝑖𝑃_𝑖+𝑐_𝑖subjectto^𝑛_𝑖=1𝑃_𝑖𝑃=𝐷_𝑖min≤𝑃_𝑖≤𝑃_𝑖^max;𝑖=1,…,𝑛, (2.1) 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, 𝑃_𝑖^min: minimum power output limit for generating unit 𝑖, and 𝑃_𝑖^max: 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 𝐶_𝑖𝑃_𝑖=𝑎_𝑖𝑃²_𝑖+𝑏_𝑖𝑃_𝑖+𝑐_𝑖+𝑒_𝑖sin𝑓_𝑖𝑃_𝑖^min−𝑃_𝑖. (2.2)

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 1Minimize2𝑥^𝑇𝑄𝑥+𝑐^𝑇𝑥subjectto𝐴𝑥=𝑏𝑙≤𝑥≤𝑢, (3.1) where 𝐴∈ℜ^𝑚×𝑛 so that rank(𝐴)=𝑛,  𝑏∈ℜ^𝑚,𝑥∈ℜ^𝑛,𝑐∈ℜ^𝑛,𝑢∈ℜ^𝑛,𝑙∈ℜ^𝑛, and𝑄∈ℜ^𝑛×𝑛 is a diagonal matrix. Problem (3.1) is equivalent to (3.2), where z∈ℜ^𝑛 is a slack variable and 𝑟∈ℜ^𝑛 is a surplus variable:1Minimize2𝑥^𝑇𝑄𝑥+𝑐^𝑇𝑥subjectto𝐴𝑥=𝑏𝑥+𝑧=𝑢;𝑧≥0𝑥−𝑟=𝑙;𝑟≥0. (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:Minimize𝐹_𝜇()=1𝑥,𝑟,𝑧2𝑥^𝑇𝑄𝑥+𝑐^𝑇𝑥−𝜇^𝑛_𝑗=1ln𝑟_𝑗−𝜇^𝑛_𝑗=1ln𝑧_𝑗subjectto𝐴𝑥=𝑏𝑥+𝑧=𝑢𝑥−𝑟=𝑙. (3.3)

The Lagrangian function related to (3.3) is expressed by𝐿_𝜇=1𝑥,𝑟,𝑧,𝑤,𝑠,𝑦2𝑥^𝑇𝑄𝑥+𝑐^𝑇𝑥+𝑤^𝑇(𝑏−𝐴𝑥)+𝑠^𝑇(𝑙−𝑥+𝑟)+𝑦^𝑇(𝑥+𝑧−𝑢)−𝜇^𝑛_𝑗=1ln𝑠_𝑗−𝜇^𝑛_𝑗=1ln𝑦_𝑗, (3.4) 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𝜕𝐿_𝜇𝜕𝑥=0⟺−𝑄𝑥−𝑐+𝐴^𝑇𝑤+𝑠−𝑦=0,𝜕𝐿_𝜇𝜕𝑤=0⟺𝐴𝑥=𝑏,𝜕𝐿_𝜇𝜕𝑦=0⟺𝑥+𝑧=𝑢,𝜕𝐿_𝜇𝜕𝑠=0⟺𝑥−𝑟=𝑙,𝜕𝐿_𝜇𝜕𝑧=0⟺𝑍𝑌𝑒−𝜇𝑒=0,𝜕𝐿_𝜇𝜕𝑟=0⟺𝑅𝑆𝑒−𝜇𝑒=0, (3.5) where 𝑅,𝑍,𝑆 and 𝑌 are diagonal matrices whose diagonal elements are 𝑟_𝑖,𝑧_𝑖,𝑠_𝑖, and 𝑦_𝑖;𝑖=1,…,𝑛, respectively; 𝑒=(1,…,1)^𝑇; 𝜇 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:Ω⁰=𝑥,𝑤,𝑧,𝑟,𝑦,𝑠⟂−𝑄𝑥+𝐴^𝑇.𝑤+𝑠−𝑦=𝑐,𝐴𝑥=𝑏,𝑥+𝑧=𝑢,𝑥−𝑟=𝑙,𝑧,𝑟,𝑦,𝑠>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.

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 ℎ^𝑘+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 𝑘+1, the new point ℎ^𝑘+1 is defined by the following equation:ℎ^𝑘+1=𝑥^𝑘+1𝑤^𝑘+1𝑧^𝑘+1𝑟^𝑘+1𝑦^𝑘+1𝑠^𝑘+1=𝑥^𝑘+𝑑^𝑘_𝑥𝑤^𝑘+𝑑^𝑘_𝑤𝑧^𝑘+𝑑^𝑘_𝑧𝑟^𝑘+𝑑^𝑘_𝑟𝑦^𝑘+𝑑^𝑘_𝑦𝑠^𝑘+𝑑^𝑘_𝑠. (3.7)

Therefore, it is necessary to determine the direction of the movement 𝑑^𝑘 to obtain the new point ℎ^𝑘+1. 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):𝐽ℎ^(𝑘)𝑑^(𝑘)=−𝐹ℎ^(𝑘),𝐴00000−𝑄𝐴^𝑇𝑑00−𝐼𝐼𝐼0𝐼000𝐼00−𝐼0000𝑌0𝑍0000𝑆0𝑅^𝑘𝑥𝑑^𝑘_𝑤𝑑^𝑘_𝑧𝑑^𝑘_𝑟𝑑^𝑘_𝑦𝑑^𝑘_𝑠=𝑡^𝑘𝑔^𝑘𝑓^𝑘𝑜^𝑘𝑞^𝑘𝑣^𝑘, (3.8)(3.9) where the residuals of (3.9) for the predictor step are expressed as in𝑡^𝑘=𝑏−𝐴𝑥^𝑘;𝑔^𝑘=𝑄𝑥^𝑘+𝑐−𝐴^𝑇𝑤^𝑘−𝑠^𝑘+𝑦^𝑘;𝑓^𝑘=𝑢−𝑥^𝑘−𝑧^𝑘;𝑜^𝑘=𝑙−𝑥^𝑘+𝑟^𝑘;𝑞^𝑘=𝜇_𝑘𝑒−𝑍_𝑘𝑌_𝑘𝑒;𝑣^𝑘=𝜇_𝑘𝑒−𝑅_𝑘𝑆_𝑘𝑒;𝑒=(1,1,…,1)^𝑇. (3.10)

The directions 𝑑^𝑘_𝑥,𝑑^𝑘_𝑤,𝑑^𝑘_𝑧,𝑑^𝑘_𝑟,𝑑^𝑘_𝑦,𝑑^𝑘_𝑠 to be determined in the predictor step use the residuals defined in (3.10), resulting in (3.11)–(3.16), as follows:𝐴𝑑^𝑘_𝑥=𝑡^𝑘,−𝑄𝑑^𝑘_𝑥+𝐴^𝑇𝑑^𝑘_𝑤+𝑑^𝑘_𝑠−𝑑^𝑘_𝑦=𝑔^𝑘𝑑_𝑘𝑥,+𝑑^𝑘_𝑧=𝑓^𝑘𝑑_𝑘𝑥,−𝑑^𝑘_𝑟=𝑜^𝑘𝑌_𝑘,𝑑^𝑘_𝑧+𝑍_𝑘𝑑^𝑘_𝑦=𝑞^𝑘𝑆_𝑘,𝑑^𝑘_𝑟+𝑅_𝑘𝑑^𝑘_𝑠=𝑣^𝑘. (3.11)(3.12)(3.13)(3.14)(3.15)(3.16)

From (3.13) and (3.14), we obtain (3.17) and (3.18), respectively:𝑑^𝑘_𝑧=−𝑑^𝑘_𝑥+𝑓^𝑘,𝑑^𝑘_𝑟=𝑑^𝑘_𝑥−𝑜^𝑘. (3.17)(3.18)

Isolating 𝑑^𝑘_𝑦 and 𝑑^𝑘_𝑠 in (3.15) and (3.16), respectively, leads to (3.19) and (3.20), as follows:𝑑^𝑘_𝑦=𝑍_𝑘⁻¹𝑞^𝑘−𝑌_𝑘𝑑^𝑘_𝑧,𝑑^𝑘_𝑠=𝑅_𝑘⁻¹𝑣^𝑘−𝑆_𝑘𝑑^𝑘_𝑟. (3.19)(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):𝜃=𝑅_𝑘⁻¹𝑆_𝑘+𝑍_𝑘⁻¹𝑌_𝑘+𝑄⁻¹,𝑑^𝑘_𝑤=𝐴𝜃𝐴^𝑇−1𝐴𝜃𝑔^𝑘+𝑝^𝑘+𝑡^𝑘𝑑_𝑘𝑥,=𝜃𝐴^𝑇𝑑^𝑘_𝑤−𝑔^𝑘−𝑝^𝑘, (3.21)(3.22)(3.23) where𝑝^𝑘=𝑅_𝑘⁻¹𝑆_𝑘𝑜^𝑘−𝑣^𝑘+𝑍_𝑘⁻¹𝑞^𝑘−𝑌_𝑘𝑓^𝑘. (3.24)

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:𝐽ℎ^𝑘𝑑^𝑘=−𝐹ℎ^𝑘, (3.25) 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𝐴00000−𝑄𝐴^𝑇𝑑00−𝐼𝐼𝐼0𝐼000𝐼00−𝐼0000𝑌0𝑍0000𝑆0𝑅^𝑘𝑥𝑑^𝑘_𝑤𝑑^𝑘_𝑧𝑑^𝑘_𝑟𝑑^𝑘_𝑦𝑑^𝑘_𝑠=𝑡^𝑘̃𝑔^𝑘𝑓^𝑘𝑜^𝑘̃𝑞^𝑘𝑣^𝑘, (3.26) where𝑡^𝑘̃=𝑏−𝐴𝑥^𝑘;𝑔^𝑘=𝑄𝑥^𝑘+𝑐−𝐴^𝑇𝑤^𝑘−𝑠^𝑘+𝑦^𝑘;𝑓^𝑘=𝑢−𝑥^𝑘−𝑧^𝑘;𝑜^𝑘̃=𝑙−𝑥^𝑘+𝑟^𝑘;𝑞^𝑘=𝜇_𝑘𝑒−𝑍_𝑘𝑌_𝑘𝑒−𝐷^𝑘_𝑧𝐷^𝑘_𝑦𝑒;𝑣^𝑘=𝜇_𝑘𝑒−𝑅_𝑘𝑆_𝑘𝑒−𝐷^𝑘_𝑥𝐷^𝑘_𝑠𝑒, (3.27)𝐷^𝑘_𝑥, 𝐷^𝑘_𝑧, 𝐷^𝑘_𝑦 and 𝐷^𝑘_𝑠 are diagonal matrices whose diagonal components are (𝑑^𝑘_𝑥)_𝑖, (𝑑^𝑘_𝑧)_𝑖, (𝑑^𝑘_𝑦)_𝑖 and (𝑑^𝑘_𝑠)_𝑖, 𝑖=1,…,𝑛, 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:𝑑^𝑘_𝑤=𝐴𝜃𝐴^𝑇−1𝐴𝜃𝑔^𝑘+𝑝^𝑘+𝑡^𝑘̃𝑑_𝑘𝑥,=𝜃𝐴^𝑇𝑑^𝑘_𝑤−𝑔^𝑘−𝑝^𝑘𝑑_𝑘𝑧,=−𝑑^𝑘_𝑥+𝑓^𝑘𝑑_𝑘𝑟,=𝑑^𝑘_𝑥+𝑜^𝑘̃𝑑_𝑘𝑦,=𝑍_𝑘⁻¹𝑞^𝑘−𝑌_𝑘𝑑^𝑘_𝑧𝑑_𝑘𝑠,=𝑅_𝑘⁻¹𝑣^𝑘−𝑆_𝑘𝑑^𝑘_𝑟, (3.28)(3.29)(3.30)(3.31)(3.32)(3.33) where𝑝^𝑘=𝑅_𝑘⁻¹𝑆_𝑘𝑜^𝑘̃−𝑣^𝑘+𝑍_𝑘⁻¹𝑞^𝑘−𝑌_𝑘𝑓^𝑘, (3.34) 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 (𝑥^𝑘+1,𝑤^𝑘+1,𝑧^𝑘+1,𝑟^𝑘+1,𝑦^𝑘+1,𝑠^𝑘+1), while ensuring that 𝑠^𝑘+1>0, 𝑧^𝑘+1>0, 𝑦^𝑘+1>0 and 𝑟^𝑘+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.

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:𝑥^𝑘+1=𝑥^𝑘+𝛼_𝑃𝑘𝑑^𝑘_𝑥,𝑟^𝑘+1=𝑟^𝑘+𝛼_𝑃𝑘𝑑^𝑘_𝑟𝑧_𝑘+1,=𝑧^𝑘+𝛼_𝑃𝑘𝑑^𝑘_𝑧𝑤_𝑘+1,=𝑤^𝑘+𝛼_𝐷𝑘𝑑^𝑘_𝑤𝑦_𝑘+1,=𝑦^𝑘+𝛼_𝐷𝑘𝑑^𝑘_𝑦𝑠_𝑘+1,=𝑠^𝑘+𝛼_𝐷𝑘𝑑^𝑘_𝑠. (3.35)(3.36)(3.37)(3.38)(3.39)(3.40)

The step size 𝛼_𝑃𝑘 for the primal variables is calculated through𝛼_𝑃𝑘=Min𝛼_𝑃,𝛼_𝑄,𝛼_LS, (3.41) where, for 0<𝛼<1, the step sizes 𝛼_𝑃,𝛼_𝑄,𝛼_LS are determined as follows.(i)The step size for primal variables 𝛼_𝑃 is obtained without violating the nonnegativity requirements of the primal variables:𝛼_𝑃=Min−𝛼𝑧^𝑘_𝑖𝑑^𝑘_𝑧𝑖,−𝛼𝑟^𝑘_𝑖𝑑^𝑘_𝑟𝑖suchthat𝑑^𝑘_𝑧𝑖,𝑑^𝑘_𝑟𝑖<0. (3.42)(ii)The step size 𝛼_𝑄 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: 𝛼_LS=Min_𝛼𝑘𝐹_𝜇𝑥^𝑘+1𝛼_𝑘=Min_𝛼𝑘𝐹_𝜇𝑥^𝑘+𝛼_𝑘𝑑^𝑘𝑥, (3.44) where 𝐹_𝜇 is defined in (3.3).

For notation simplicity, in the summary of the Fibonacci search, the following identities are defined: 𝑥^𝑘+1(𝛼_LS)=(𝑥^𝑘+1(𝛼_LS),𝑟^𝑘+1(𝛼_LS),𝑧^𝑘+1(𝛼_LS))^𝑇=(𝑥^𝑘+𝛼_LS𝑑^𝑘_𝑥,𝑟^𝑘+𝛼_LS𝑑^𝑘_𝑟,𝑧^𝑘+𝛼_LS𝑑^𝑘_𝑧)^𝑇, 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 𝑥^𝑘+1 in direction 𝑑𝑥^𝑘 using the function 𝐹_𝜇 defined in (3.3). The Fibonacci method calculates 𝛼_LS so that the minimization of 𝐹_𝜇 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 𝛼_𝐷𝑘 for the dual variables (3.38)–(3.40) is calculated through (3.45) for 0<𝛼<1:𝛼_𝐷𝑘=Min1,−𝛼𝑠^𝑘_𝑖𝑑^𝑘_𝑠𝑖,−𝛼𝑦^𝑘_𝑖𝑑^𝑘_𝑦𝑖suchthat𝑑^𝑘_𝑠𝑖;𝑑^𝑘_𝑦𝑖<0. (3.45)

3.4.2. Corrector Step Size

In the corrector step, the calculation of step size is defined analogously, considering the variables from the following 𝑥^𝑘+1=𝑥^𝑘+𝛽_𝑃𝑘𝑑^𝑘_𝑥,𝑟^𝑘+1=𝑟^𝑘+𝛽_𝑃𝑘𝑑^𝑘_𝑟𝑧_𝑘+1,=𝑧^𝑘+𝛽_𝑃𝑘𝑑^𝑘_𝑧𝑤_𝑘+1,=𝑤^𝑘+𝛽_𝐷𝑘𝑑^𝑘_𝑤𝑦_𝑘+1,=𝑦^𝑘+𝛽_𝐷𝑘𝑑^𝑘_𝑦𝑠_𝑘+1,=𝑠^𝑘+𝛽_𝐷𝑘𝑑^𝑘_𝑦, (3.46)(3.47)(3.48)(3.49)(3.50)(3.51) where𝛽_𝑃𝑘=Min𝛽_𝑃,𝛽_𝑄,𝛽_LS, (3.52) while the step sizes 𝛽_𝑃,𝛽_𝑄,𝛽_LS are determined as follows.(i)The step size for primal variables 𝛽_𝑃 is obtained without violating the nonnegativity requirements of the primal variables:𝛽_𝑃=Min1,−𝛼𝑧^𝑘_𝑖𝑑^𝑘_𝑧𝑖,−𝛼𝑟^𝑘_𝑖𝑑^𝑘_𝑟𝑖suchthat𝑑^𝑘_𝑧𝑖,𝑑^𝑘_𝑟𝑖<0. (3.53)(ii)The step size 𝛽_𝑄 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:𝛽_LS=𝛽_𝑘, (3.55) so that𝐹_𝜇𝑥^𝑘+1𝛽_𝑘≤𝐹_𝜇𝑥^𝑘+𝛽_𝑘∇𝐹_𝜇𝑥^𝑘𝑇𝑑^𝑘𝑥, (3.56) where 𝑥^𝑘+1(𝛽_𝑘)=(𝑥^𝑘+1(𝛽_𝑘),𝑟^𝑘+1(𝛽_𝑘),𝑧^𝑘+1(𝛽_𝑘))^𝑇=(𝑥^𝑘+𝛽_𝑘𝑑^𝑘_𝑥,𝑟^𝑘+𝛽_𝑘𝑑^𝑘_𝑟,𝑧^𝑘+𝛽_𝑘𝑑^𝑘_𝑧)^𝑇; 𝑑𝑥^𝑘=(𝑑^𝑘_𝑥,𝑑^𝑘_𝑟,𝑑^𝑘_𝑧); 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 𝑥^𝑘+1, in direction 𝑑𝑥^𝑘, so that the function 𝐹_𝜇 defined in (3.3) decreases. This search calculates 𝛽_LS, 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 𝛽₀=1. If (3.44) is not satisfied, 𝛽_𝑘 is updated using the following sequence:𝛽_𝑘+1=𝛽_𝑘𝛿,𝑘=0,1,…,𝑛;with𝛿>1,usuallysetto𝛿=2. (3.57)

The following equation is used for the analysis of (3.56): ∇𝐹_𝜇𝑥^𝑘=∇𝐹_𝜇𝑥^𝑘,∇𝐹_𝜇𝑟^𝑘,∇𝐹_𝜇𝑧^𝑘𝑇=𝑄𝑥^𝑘+𝑐,−𝜇_𝑘𝑅_𝑘⁻¹𝑒,−𝜇_𝑘𝑍_𝑘⁻¹𝑒^𝑇 (3.58) 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:𝛽_𝐷𝑘=Min1,−𝛼𝑠^𝑘_𝑖𝑑^𝑘_𝑠𝑖,−𝛼𝑦^𝑘_𝑖𝑑^𝑘_𝑦𝑖suchthat𝑑^𝑘_𝑠𝑖;𝑑^𝑘_𝑦𝑖<0. (3.59)

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:𝑡^𝑘‖𝑏‖+1≤𝜀₁,𝑡̃predictorstep_𝑘‖𝑏‖+1≤𝜀₁correctorstep; (3.60)(3.61)

(ii)dual feasibility:𝑢^𝑘𝑄𝑥^𝑘+𝑐+1≤𝜀₂predictorstep,𝑢^𝑘𝑄𝑥^𝑘+𝑐+1≤𝜀₂correctorstep; (3.62)(3.63)

(iii)complementary slackness:𝑣^𝑘≤𝜀₃,𝑞^𝑘≤𝜀₃predictorstep𝑣^𝑘≤𝜀₃,𝑞^𝑘≤𝜀₃correctorstep, (3.64)(3.65) 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𝜇¹_𝑘=𝑟^𝑘𝑇𝑠^𝑘𝑛;𝜇²_𝑘=𝑧^𝑘𝑇𝑦^𝑘𝑛, (3.66) so that 𝜇_𝑘 is calculated using the following equation𝜇_𝑘=Min𝜎𝜇¹_𝑘,𝜎𝜇²_𝑘foraconstant0<𝜎<1, (3.67) 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 𝑘=0. 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 𝑥^𝑘+1; 𝑤^𝑘+1; 𝑠^𝑘+1; 𝑦^𝑘+1; 𝑧^𝑘+1;𝑟^𝑘+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: 𝑡^𝑘̃; 𝑔^𝑘; 𝑓^𝑘; 𝑜^𝑘̃; 𝑞^𝑘; 𝑣^𝑘 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 𝑥^𝑘+1; 𝑤^𝑘+1; 𝑠^𝑘+1; 𝑦^𝑘+1; 𝑧^𝑘+1; 𝑟^𝑘+1 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.

Step 12. Adjust 𝑘=𝑘+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.

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 𝑃_𝑖^min,𝑃_𝑖^max represent minimum and maximum power output capabilities, respectively, for generating unit 𝑖.

The following values are adopted to initialise the method:𝑥⁰=(450,100,300);𝑤⁰=(0,0,0);𝑦⁰=(0,0,0). (4.1)

The parameters related to the system’s total demand are set at 𝐷=850MW and the active power losses are neglected. The values adopted for 𝜀₁,𝜀₂, and 𝜀₃ are 10⁻⁸.

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 𝑃_𝑖^min,𝑃_𝑖^max represent minimum and maximum power output capabilities, respectively, for the generating unit 𝑖.

The following values are adopted to initialise the method:𝑥⁰=(20,30,75,75,145,155);𝑤⁰=(0,0,0,0,0,0);𝑦⁰=(1,1,1,1,1,1). (4.2)

The parameters related to the system’s total demand are set at 𝐷=500MW and the active power losses are neglected. The values adopted for 𝜀₁,𝜀₂, and 𝜀₃ are 10⁻⁸.

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:𝑥⁰=();𝑤660,320,330,150,140,155,165,150,140,70,80,75,85₀=();𝑦0,0,0,0,0,0,0,0,0,0,0,0,0₀=(30,30,30,30,30,30,30,30,30,30,30,30,30). (4.3)

The parameters related to the system’s total demand are set at 𝐷=2520MW and the active power losses are neglected. The values adopted for 𝜀₁,𝜀₂, and 𝜀₃ are 10⁻⁸.

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:Minimize𝐶_𝑇subjectto^𝑛_𝑖=1𝑃_𝑖𝑃=𝐷_𝑖∈𝑆, (4.4) where 𝑆={𝑃_𝑖∈ℜ⟂𝑃_𝑖−𝑃_𝑖^min≥0and𝑃_𝑖^max−𝑃_𝑖≥0,𝑖=1…,𝑛}, and 𝐶_𝑇 = ∑^𝑛_𝑖=1(𝐶_𝑖(𝑃_𝑖) = 𝑎_𝑖𝑃²_𝑖+𝑏_𝑖𝑃_𝑖+𝑐_𝑖).

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:𝐿𝑃_𝑖,𝑤,𝑠_𝑖,𝑦_𝑖=𝐶_𝑇𝑃_𝑖+𝑤𝐷−^𝑛_𝑖=1𝑃_𝑖+^𝑛_𝑖=1𝑠_𝑖𝑃_𝑖^min−𝑃_𝑖+^𝑛_𝑖=1𝑦_𝑖𝑃_𝑖−𝑃_𝑖^max, (4.5) where 𝑤∈ℜ,  𝑠_𝑖∈ℜ, and 𝑦_𝑖∈ℜ,  𝑖=1,…,𝑛, are the associated Lagrange multipliers. The KKT conditions associated with the Lagrangian function (4.5) are expressed through𝐷−^𝑛_𝑖=1𝑃_𝑖𝑃=0,_𝑖−𝑃_𝑖^min≥0;𝑠_𝑖𝑃_𝑖−𝑃_𝑖^min𝑃=0;𝑖=1,…,𝑛,_𝑖max−𝑃_𝑖≥0;𝑦_𝑖𝑃_𝑖^max−𝑃_𝑖𝜆=0;𝑖=1,…,𝑛,_𝑖𝑃_𝑖−𝑤−𝑠_𝑖+𝑦_𝑖=0;𝑖=1,…,𝑛, (4.6)(4.7)(4.8)(4.9) where𝜆_𝑖𝑃_𝑖=𝜕𝐶_𝑖𝑃_𝑖𝜕𝑃_𝑖=2𝑎_𝑖𝑃_𝑖+𝑏_𝑖;𝑖=1,…,𝑛. (4.10)

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𝑤=𝜆_𝑖𝑃_𝑖−𝑠_𝑖+𝑦_𝑖;𝑖=1,…,𝑛. (4.11)

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𝑠^∗_𝑖=0, 𝑦^∗_𝑖=0; 𝑖=1,…,𝑛. Thus, from (4.11), one reaches𝑤=𝜆^∗_𝑖𝑃^∗_𝑖;𝑖=1,…,𝑛, (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.

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Ω_max be defined with the indices of the generating units that have been optimally dispatched in their maximum power output. In this situation, 𝑠^∗_𝑖=0 for 𝑖=1,…,𝑛, 𝑦^∗_𝑗>0 for 𝑗∈Ω_max, and 𝑦^∗_𝑗=0for 𝑗∉Ω_max. Starting from (4.11), one reaches𝑤=𝜆^∗_𝑗𝑃^∗_𝑗+𝑦^∗_𝑗;𝑗∈Ω_max,𝑤=𝜆^∗_𝑗𝑃^∗_𝑗;𝑗∉Ω_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 𝑦^∗_𝑗>0;𝑗∈Ω_max.

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 Ω_min be defined with the indices of the generating units that have been optimally dispatched in their minimum power output. In this situation, 𝑦^∗_𝑖=0 for 𝑖=1,…,𝑛, 𝑠^∗_𝑗>0 for 𝑗∈Ω_min, and 𝑠^∗_𝑗=0 for 𝑗∉Ω_min. Thus, using (4.11), one reaches𝑤=𝜆^∗_𝑗𝑃^∗_𝑗−𝑠^∗_𝑗;𝑗∈Ω_min,𝑤=𝜆^∗_𝑗𝑃^∗_𝑗;𝑗∉Ω_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 𝑠^∗_𝑗>0;𝑗∈Ω_min.

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, 𝑃^∗_𝑖,𝑖=1,…,𝑛, 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 𝑦^∗_𝑖;𝑖=1,…,𝑛 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, 𝑃^∗_𝑖,𝑖=1,…,𝑛, 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 𝑦^∗_𝑖;𝑖=1,…,𝑛 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, 𝑃^∗_𝑖,𝑖=1,…,𝑛, 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𝑦^∗_𝑖; 𝑖=1,…,𝑛 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 𝑃^∗_𝑖,𝑖=1,…,𝑛 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).