# Predictor-Corrector Primal-Dual Interior Point Method for Solving Economic

Authors: Antonio Roberto Balbo; Márcio Augusto da Silva Souza; Edméa Cássia Baptista; Leonardo Nepomuceno
Journal: Mathematical Problems in Engineering
 (2012)

---

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


---

## Body

## 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: Section2 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(2.1)MinimizeCT=∑i=1n(Ci(Pi)=aiPi2+biPi+ci)subjectto∑i=1nPi=DPimin⁡≤Pi≤Pimax⁡;i=1,…,n,
where CT: total cost function for generating units, n: number of generating units, Ci(Pi): cost of generation unit i (without considering the valve-point effect), ai,bi,ci: coefficient of the cost function for generating unit i, Pi: power output of generating unit i, D: total power demand, Pimin: minimum power output limit for generating unit i, and Pimax: maximum power output limit for generating unit i.The objective function in (2.1) may also represent the so-called valve-point effects, as in [28], which are associated with the opening of pressure valves at some specific operating points. In such cases, Ci(Pi) is mathematically described as in (2.2)Ci(Pi)=aiPi2+biPi+ci+|eisin⁡(fi(Pimin⁡-Pi))|.The cost functionCi(Pi) 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.


## 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(2.1)MinimizeCT=∑i=1n(Ci(Pi)=aiPi2+biPi+ci)subjectto∑i=1nPi=DPimin⁡≤Pi≤Pimax⁡;i=1,…,n,
where CT: total cost function for generating units, n: number of generating units, Ci(Pi): cost of generation unit i (without considering the valve-point effect), ai,bi,ci: coefficient of the cost function for generating unit i, Pi: power output of generating unit i, D: total power demand, Pimin: minimum power output limit for generating unit i, and Pimax: maximum power output limit for generating unit i.The objective function in (2.1) may also represent the so-called valve-point effects, as in [28], which are associated with the opening of pressure valves at some specific operating points. In such cases, Ci(Pi) is mathematically described as in (2.2)Ci(Pi)=aiPi2+biPi+ci+|eisin⁡(fi(Pimin⁡-Pi))|.The cost functionCi(Pi) 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
(3.1)Minimize12xTQx+cTxsubject toAx=bl≤x≤u,
where A∈ℜm×n so that rank(A)=n,  b∈ℜm,x∈ℜn,c∈ℜn,u∈ℜn,l∈ℜn, andQ∈ℜn×n is a diagonal matrix. Problem (3.1) is equivalent to (3.2), where z∈ℜn is a slack variable and r∈ℜn is a surplus variable:(3.2)Minimize12xTQx+cTxsubject toAx=bx+z=u;z≥0x-r=l;r≥0.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:(3.3)MinimizeFμ(x,r,z)=12xTQx+cTx-μ∑j=1nln⁡rj-μ∑j=1nln⁡zjsubjecttoAx=bx+z=ux-r=l.The Lagrangian function related to (3.3) is expressed by(3.4)Lμ(x,r,z,w,s,y)=12xTQx+cTx+wT(b-Ax)+sT(l-x+r)+yT(x+z-u)-μ∑j=1nln⁡⁡sj-μ∑j=1nln⁡yj,
where the dual variables associated with the three equality constraints in (3.3) are, respectively, w∈ℜm,y∈ℜn,s∈ℜn. The optimal Karush-Kuhn-Tucker (KKT) conditions for problem (3.1) are depicted in(3.5)∂Lμ∂x=0⟺-Qx-c+ATw+s-y=0,∂Lμ∂w=0⟺Ax=b,∂Lμ∂y=0⟺x+z=u,∂Lμ∂s=0⟺x-r=l,∂Lμ∂z=0⟺ZYe-μe=0,∂Lμ∂r=0⟺RSe-μe=0,
where R,Z,S and Y are diagonal matrices whose diagonal elements are ri,zi,si, and yi;i=1,…,n, respectively; e=(1,…,1)T; μ is a dual metric or adjustment parameter for the curve defined by the central trajectory (path-following parameter). The set Ω0 given in (3.6) is defined to simplify to notation. This set describes interior points for problem (3.2) and its corresponding dual problem:(3.6)Ω0={(x,w,z,r,y,s)⊥-Qx+ATw+s-y=c,Ax=b,x+z=u,x-r=l,(z,r,y,s)>0AT(x,w,z,r,y,s)}.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 Section3.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 iterationk, a point hk satisfies the primal and dual feasibility conditions expressed by (3.5). In this case, the definition of the new point hk+1 depends solely on the calculation of a search direction and on a step size in such a direction. Disregarding the step size in the iteration k+1, the new point hk+1 is defined by the following equation:(3.7)hk+1=(xk+1wk+1zk+1rk+1yk+1sk+1)=(xk+dxkwk+dwkzk+dzkrk+drkyk+dyksk+dsk).Therefore, it is necessary to determine the direction of the movementdk to obtain the new point hk+1. Following the steps in the Newton method applied to the nonlinear system (3.5), dk can be obtained by solving for the system (3.8), which is equivalent to  (3.9):(3.8)J(h(k))d(k)=-F(h(k)),(3.9)(A00000-QAT00-III0I000I00-I0000Y0Z0000S0R)(dxkdwkdzkdrkdykdsk)=(tkgkfkokqkvk),
where the residuals of (3.9) for the predictor step are expressed as in(3.10)tk=b-Axk;gk=Qxk+c-ATwk-sk+yk;fk=u-xk-zk;ok=l-xk+rk;qk=μke-ZkYke;vk=μke-RkSke;e=(1,1,…,1)T.The directionsdxk,dwk,dzk,drk,dyk,dsk to be determined in the predictor step use the residuals defined in (3.10), resulting in (3.11)–(3.16), as follows:(3.11)Adxk=tk,(3.12)-Qdxk+ATdwk+dsk-dyk=gk,(3.13)dxk+dzk=fk,(3.14)dxk-drk=ok,(3.15)Ykdzk+Zkdyk=qk,(3.16)Skdrk+Rkdsk=vk.From (3.13) and (3.14), we obtain (3.17) and (3.18), respectively:(3.17)dzk=-dxk+fk,(3.18)drk=dxk-ok.Isolatingdyk and dsk in (3.15) and (3.16), respectively, leads to (3.19) and (3.20), as follows:(3.19)dyk=Zk-1(qk-Ykdzk),(3.20)dsk=Rk-1(vk-Skdrk).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):(3.21)θ=(Rk-1Sk+Zk-1Yk+Q)-1,(3.22)dwk=(AθAT)-1[Aθ(gk+pk)+tk],(3.23)dxk=θ(ATdwk-gk-pk),
where(3.24)pk=Rk-1(Skok-vk)+Zk-1(qk-Ykfk).Note that, after calculating (3.23), the remaining components of the direction vector dzk,drk,dyk, and dsk in (3.17), (3.18), (3.19), and (3.20), respectively, are easily calculated. Since matrix AθAT is symmetrical and positive definite in (3.22) (considering that Q is symmetrical and positive definite), dwk can be determined by using the Cholesky decomposition.


### 3.3. Search Directions—Corrector Step

Analogously to the procedure carried out in Section3.2, this section describes the calculation of the search direction for the corrector stepd̃k, which is obtained by solving the linear system (3.25), as follows:(3.25)J(hk)d̃k=-F̃(hk),
where F̃(hk) is obtained by considering second-order approximations in the residuals (3.10) (calculated in the predictor step), so that (3.25) is equivalent to(3.26)(A00000-QAT00-III0I000I00-I0000Y0Z0000S0R)(d̃xkd̃wkd̃zkd̃rkd̃ykd̃sk)=(t̃kg̃kf̃kõkq̃kṽk),
where(3.27)t̃k=b-Axk;g̃k=Qxk+c-ATwk-sk+yk;f̃k=u-xk-zk;õk=l-xk+rk;q̃k=μke-ZkYke-DzkDyke;ṽk=μke-RkSke-DxkDske,Dxk, Dzk, Dyk and Dsk are diagonal matrices whose diagonal components are (dxk)i, (dzk)i, (dyk)i and (dsk)i, i=1,…,n, respectively.The calculation of the residualstk,gk,fk,ok,qk,vk, described in (3.10), and of t̃k,g̃k,f̃k,õk,q̃k,ṽk, described in (3.27), basically distinguishes the predictor and corrector steps in the proposed method. It is important to note that the corrector step procedure uses direction values dxk, dzk, dyk, and dsk, which have already been calculated in the predictor step, to redefine residuals q̃k and ṽk in (3.27). Using (3.25) and (3.26) and following the same steps taken to determine the directions of the predictor step, as seen in Section 3.2, the components of the direction vector d̃kcan be calculated using the following:(3.28)d̃wk=(AθAT)-1[Aθ(gk+p̃k)+t̃k],(3.29)d̃xk=θ(ATd̃wk-gk-p̃k),(3.30)d̃zk=-d̃xk+f̃k,(3.31)d̃rk=d̃xk+õk,(3.32)d̃yk=Zk-1(q̃k-Ykd̃zk),(3.33)d̃sk=Rk-1(ṽk-Skd̃rk),
where(3.34)p̃k=Rk-1(Skõk-ṽk)+Zk-1(q̃k-Ykf̃k),
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(xk+1,wk+1,zk+1,rk+1,yk+1,sk+1), while ensuring that sk+1>0, zk+1>0, yk+1>0 and rk+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:(3.35)xk+1=xk+αPkdxk,(3.36)rk+1=rk+αPkdrk,(3.37)zk+1=zk+αPkdzk,(3.38)wk+1=wk+αDkdwk,(3.39)yk+1=yk+αDkdyk,(3.40)sk+1=sk+αDkdsk.The step sizeαPk for the primal variables is calculated through(3.41)αPk=Min⁡{αP,αQ,αLS},
where, for 0<α<1, the step sizes αP,αQ,αLS are determined as follows.(i)
The step size for primal variablesαP is obtained without violating the nonnegativity requirements of the primal variables:(3.42)αP=Min⁡{-αzikdzik,-αrikdriksuch thatdzik,drik<0}.(ii)
The step sizeαQ is the maximum step size possible without increasing the objective value:(3.43)αQ=-[(dxk)t(Qxk+c)][(dxk)tQdxk].(iii)
The step sizeαLS is determined as in (3.44) from the Fibonacci line search strategy, which is briefly summarized as follows:
(3.44)αLS=Min⁡αk⁡{Fμ(x̂k+1(αk))}=Min⁡αk⁡{Fμ(x̂k+αkdx̂k)},
where Fμ is defined in (3.3).For notation simplicity, in the summary of the Fibonacci search, the following identities are defined:x̂k+1(αLS)=(xk+1(αLS),rk+1(αLS),zk+1(αLS))T=(xk+αLSdxk,rk+αLSdrk,zk+αLSdzk)T, and dx̂k=(dxk,drk,dzk); so that dxkis defined in (3.23), drkis defined in (3.18), and dzkis defined in (3.17). Only primal variables are considered in the Fibonacci search algorithm used by the IPLS method. Starting at a point x̂k=(xk,rk,zk)T, the algorithm searches a new point x̂k+1 in direction dx̂k using the function Fμ defined in (3.3). The Fibonacci method calculates αLS so that the minimization of Fμ is ensured. αLS is determined considering the Fibonacci sequence. The initial value of αLS is set taking into account the interval of uncertainty [0,1].The step sizeαDk for the dual variables (3.38)–(3.40) is calculated through (3.45) for 0<α<1:(3.45)αDk=Min⁡{1,-αsikdsik,-αyikdyiksuch thatdsik;dyik<0}.


#### 3.4.2. Corrector Step Size

In the corrector step, the calculation of step size is defined analogously, considering the variables from the following(3.46)xk+1=xk+βPkd̃xk,(3.47)rk+1=rk+βPkd̃rk,(3.48)zk+1=zk+βPkd̃zk,(3.49)wk+1=wk+βDkd̃wk,(3.50)yk+1=yk+βDkd̃yk,(3.51)sk+1=sk+βDkd̃yk,
where(3.52)βPk=Min⁡{βP,βQ,βLS},
while the step sizes βP,βQ,βLS are determined as follows.(i)
The step size for primal variablesβP is obtained without violating the nonnegativity requirements of the primal variables:(3.53)βP=Min⁡{1,-αzikd̃zik,-αrikd̃riksuch thatd̃zik,d̃rik<0}.(ii)
The step sizeβQ is the maximum step size possible without increasing the objective value:(3.54)βQ=-[(d̃xk)t(Qxk+c)][(d̃xk)tQdxk].(iii)
The step sizeβLS in (3.55) is determined from the Armijo line search:(3.55)βLS=βk,
so that(3.56)Fμ(x̂k+1(βk))≤Fμ(x̂k)+βk∇Fμ(x̂k)Td̃x̂k,
where x̂k+1(βk)=(xk+1(βk),rk+1(βk),zk+1(βk))T=(xk+βkd̃xk,rk+βkd̃rk,zk+βkd̃zk)T; d̃x̂k=(d̃xk,d̃rk,d̃zk); so that d̃xk,d̃rk,d̃zkare defined in (3.29)–(3.31), respectively. In the proposed IPLS method, the Armijo search is also performed in the direction defined by primal variables. Starting at a point x̂k=(xk,rk,zk), the method searches for a new point x̂k+1, in direction d̃x̂k, so that the function Fμ defined in (3.3) decreases. This search calculates βLS, thereby ensuring the reduction of Fμ. To prevent oscillations in the iterative process, the initial choice for βk should not be too high, and to prevent the process from stopping prematurely, it should not be too low. Therefore, this value is generally adjusted to β0=1. If (3.44) is not satisfied, βk is updated using the following sequence:(3.57)βk+1=βkδ,k=0,1,…,n;withδ>1,usually set toδ=2.The following equation is used for the analysis of (3.56):
(3.58)∇Fμ(x̂k)=(∇Fμ(xk),∇Fμ(rk),∇Fμ(zk))T=(Qxk+c,-μkRk-1e,-μkZk-1e)T
with μk determined as in Section 3.6.The step size for dual variablesβDk is calculated by (3.59), without violating the nonnegativity requirements of the dual variables:(3.59)βDk=Min⁡{1,-αsikd̃sik,-αyikd̃yiksuch thatd̃sik;d̃yik<0}.The general principle used here to calculateαPk and βPk 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 residualstk,uk, and also for the dual metric μk. Nevertheless, it is possible to use relative values for the metrics tk,uk, and μk 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:(3.60)‖tk‖‖b‖+1≤ε1(predictor step),(3.61)‖t̃k‖‖b‖+1≤ε1(corrector step);(ii)
dual feasibility:(3.62)‖uk‖‖Qxk+c‖+1≤ε2(predictor step),(3.63)‖ũk‖‖Qxk+c‖+1≤ε2(corrector step);(iii)
complementary slackness:(3.64)‖vk‖≤ε3,‖qk‖≤ε3(predictor step)(3.65)‖ṽk‖≤ε3,‖q̃k‖≤ε3(corrector step),
where ε1, ε2, and ε3 are sufficiently small positive numbers. Eventually, other criteria can be adopted according to the specific characteristics of each problem, as can be seen in [24, 32].


### 3.6. Update of the Barrier Parameter

According to [32], the barrier parameter is updated using an inner product that involves the primal variables rk and zk, and the dual variables sk and yk, respectively, as shown in the following equation(3.66)μk1=(rk)Tskn;μk2=(zk)Tykn,
so that μk is calculated using the following equation(3.67)μk=Min⁡{σμk1,σμk2}for a constant0<σ<1,
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).
Adjustk=0. Choose an arbitrary point:(x0;w0;z0,y0,r0,y0,s0)∈Ω0 and choose ε1, ε2, and ε3 as sufficiently small positive numbers.Step 2 (intermediate calculations—predictor).
Calculate:tk;gk;fk;ok;qk;vk using (3.10), μk using (3.67) and matrix θ using (3.21).Step 3 (finding directions of translation—predictor).
Calculate search directionsdxk, dwk, dzk, drk, dyk, and dsk for the predictor step using (3.17)–(3.24).Step 4 (computing step size—predictor).
Calculate Fibonacci step sizesαPk and αDk using (3.41)–(3.44).Step 5 (moving to a new solution—predictor).
Updatexk+1; wk+1; sk+1; yk+1; zk+1;rk+1 obtained from the predictor step, according to (3.35)–(3.44).Step 6 (checking for optimality—predictor).
If the criteria defined in (3.60), (3.62), and (3.64) are satisfied, stop. The solution is optimal. Otherwise, go on to the following step.Step 7 (intermediate calculations—corrector).
Calculate:t̃k; g̃k; f̃k; õk; q̃k; ṽk using (3.27), μk using (3.67) and matrix θ using (3.21).Step 8 (finding directions of translation—corrector).
Calculated̃xk, d̃wk, d̃zk, d̃rk, d̃yk, and d̃sk, for the corrector step using (3.28)–(3.33).Step 9 (computing step size—corrector).
Calculate the Armijo step sizeβPk and βDk using (3.52)–(3.58).Step 10 (moving to a new solution—corrector).
Updatexk+1; wk+1; sk+1; yk+1; zk+1; rk+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.
Adjustk=k+1 and return to Step 2.The predictor Steps2 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.


## 3.1. Search Directions

In this section, the search directions used in the proposed IPLS method are investigated. In Section3.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 iterationk, a point hk satisfies the primal and dual feasibility conditions expressed by (3.5). In this case, the definition of the new point hk+1 depends solely on the calculation of a search direction and on a step size in such a direction. Disregarding the step size in the iteration k+1, the new point hk+1 is defined by the following equation:(3.7)hk+1=(xk+1wk+1zk+1rk+1yk+1sk+1)=(xk+dxkwk+dwkzk+dzkrk+drkyk+dyksk+dsk).Therefore, it is necessary to determine the direction of the movementdk to obtain the new point hk+1. Following the steps in the Newton method applied to the nonlinear system (3.5), dk can be obtained by solving for the system (3.8), which is equivalent to  (3.9):(3.8)J(h(k))d(k)=-F(h(k)),(3.9)(A00000-QAT00-III0I000I00-I0000Y0Z0000S0R)(dxkdwkdzkdrkdykdsk)=(tkgkfkokqkvk),
where the residuals of (3.9) for the predictor step are expressed as in(3.10)tk=b-Axk;gk=Qxk+c-ATwk-sk+yk;fk=u-xk-zk;ok=l-xk+rk;qk=μke-ZkYke;vk=μke-RkSke;e=(1,1,…,1)T.The directionsdxk,dwk,dzk,drk,dyk,dsk to be determined in the predictor step use the residuals defined in (3.10), resulting in (3.11)–(3.16), as follows:(3.11)Adxk=tk,(3.12)-Qdxk+ATdwk+dsk-dyk=gk,(3.13)dxk+dzk=fk,(3.14)dxk-drk=ok,(3.15)Ykdzk+Zkdyk=qk,(3.16)Skdrk+Rkdsk=vk.From (3.13) and (3.14), we obtain (3.17) and (3.18), respectively:(3.17)dzk=-dxk+fk,(3.18)drk=dxk-ok.Isolatingdyk and dsk in (3.15) and (3.16), respectively, leads to (3.19) and (3.20), as follows:(3.19)dyk=Zk-1(qk-Ykdzk),(3.20)dsk=Rk-1(vk-Skdrk).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):(3.21)θ=(Rk-1Sk+Zk-1Yk+Q)-1,(3.22)dwk=(AθAT)-1[Aθ(gk+pk)+tk],(3.23)dxk=θ(ATdwk-gk-pk),
where(3.24)pk=Rk-1(Skok-vk)+Zk-1(qk-Ykfk).Note that, after calculating (3.23), the remaining components of the direction vector dzk,drk,dyk, and dsk in (3.17), (3.18), (3.19), and (3.20), respectively, are easily calculated. Since matrix AθAT is symmetrical and positive definite in (3.22) (considering that Q is symmetrical and positive definite), dwk can be determined by using the Cholesky decomposition.


## 3.3. Search Directions—Corrector Step

Analogously to the procedure carried out in Section3.2, this section describes the calculation of the search direction for the corrector stepd̃k, which is obtained by solving the linear system (3.25), as follows:(3.25)J(hk)d̃k=-F̃(hk),
where F̃(hk) is obtained by considering second-order approximations in the residuals (3.10) (calculated in the predictor step), so that (3.25) is equivalent to(3.26)(A00000-QAT00-III0I000I00-I0000Y0Z0000S0R)(d̃xkd̃wkd̃zkd̃rkd̃ykd̃sk)=(t̃kg̃kf̃kõkq̃kṽk),
where(3.27)t̃k=b-Axk;g̃k=Qxk+c-ATwk-sk+yk;f̃k=u-xk-zk;õk=l-xk+rk;q̃k=μke-ZkYke-DzkDyke;ṽk=μke-RkSke-DxkDske,Dxk, Dzk, Dyk and Dsk are diagonal matrices whose diagonal components are (dxk)i, (dzk)i, (dyk)i and (dsk)i, i=1,…,n, respectively.The calculation of the residualstk,gk,fk,ok,qk,vk, described in (3.10), and of t̃k,g̃k,f̃k,õk,q̃k,ṽk, described in (3.27), basically distinguishes the predictor and corrector steps in the proposed method. It is important to note that the corrector step procedure uses direction values dxk, dzk, dyk, and dsk, which have already been calculated in the predictor step, to redefine residuals q̃k and ṽk in (3.27). Using (3.25) and (3.26) and following the same steps taken to determine the directions of the predictor step, as seen in Section 3.2, the components of the direction vector d̃kcan be calculated using the following:(3.28)d̃wk=(AθAT)-1[Aθ(gk+p̃k)+t̃k],(3.29)d̃xk=θ(ATd̃wk-gk-p̃k),(3.30)d̃zk=-d̃xk+f̃k,(3.31)d̃rk=d̃xk+õk,(3.32)d̃yk=Zk-1(q̃k-Ykd̃zk),(3.33)d̃sk=Rk-1(ṽk-Skd̃rk),
where(3.34)p̃k=Rk-1(Skõk-ṽk)+Zk-1(q̃k-Ykf̃k),
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(xk+1,wk+1,zk+1,rk+1,yk+1,sk+1), while ensuring that sk+1>0, zk+1>0, yk+1>0 and rk+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:(3.35)xk+1=xk+αPkdxk,(3.36)rk+1=rk+αPkdrk,(3.37)zk+1=zk+αPkdzk,(3.38)wk+1=wk+αDkdwk,(3.39)yk+1=yk+αDkdyk,(3.40)sk+1=sk+αDkdsk.The step sizeαPk for the primal variables is calculated through(3.41)αPk=Min⁡{αP,αQ,αLS},
where, for 0<α<1, the step sizes αP,αQ,αLS are determined as follows.(i)
The step size for primal variablesαP is obtained without violating the nonnegativity requirements of the primal variables:(3.42)αP=Min⁡{-αzikdzik,-αrikdriksuch thatdzik,drik<0}.(ii)
The step sizeαQ is the maximum step size possible without increasing the objective value:(3.43)αQ=-[(dxk)t(Qxk+c)][(dxk)tQdxk].(iii)
The step sizeαLS is determined as in (3.44) from the Fibonacci line search strategy, which is briefly summarized as follows:
(3.44)αLS=Min⁡αk⁡{Fμ(x̂k+1(αk))}=Min⁡αk⁡{Fμ(x̂k+αkdx̂k)},
where Fμ is defined in (3.3).For notation simplicity, in the summary of the Fibonacci search, the following identities are defined:x̂k+1(αLS)=(xk+1(αLS),rk+1(αLS),zk+1(αLS))T=(xk+αLSdxk,rk+αLSdrk,zk+αLSdzk)T, and dx̂k=(dxk,drk,dzk); so that dxkis defined in (3.23), drkis defined in (3.18), and dzkis defined in (3.17). Only primal variables are considered in the Fibonacci search algorithm used by the IPLS method. Starting at a point x̂k=(xk,rk,zk)T, the algorithm searches a new point x̂k+1 in direction dx̂k using the function Fμ defined in (3.3). The Fibonacci method calculates αLS so that the minimization of Fμ is ensured. αLS is determined considering the Fibonacci sequence. The initial value of αLS is set taking into account the interval of uncertainty [0,1].The step sizeαDk for the dual variables (3.38)–(3.40) is calculated through (3.45) for 0<α<1:(3.45)αDk=Min⁡{1,-αsikdsik,-αyikdyiksuch thatdsik;dyik<0}.


### 3.4.2. Corrector Step Size

In the corrector step, the calculation of step size is defined analogously, considering the variables from the following(3.46)xk+1=xk+βPkd̃xk,(3.47)rk+1=rk+βPkd̃rk,(3.48)zk+1=zk+βPkd̃zk,(3.49)wk+1=wk+βDkd̃wk,(3.50)yk+1=yk+βDkd̃yk,(3.51)sk+1=sk+βDkd̃yk,
where(3.52)βPk=Min⁡{βP,βQ,βLS},
while the step sizes βP,βQ,βLS are determined as follows.(i)
The step size for primal variablesβP is obtained without violating the nonnegativity requirements of the primal variables:(3.53)βP=Min⁡{1,-αzikd̃zik,-αrikd̃riksuch thatd̃zik,d̃rik<0}.(ii)
The step sizeβQ is the maximum step size possible without increasing the objective value:(3.54)βQ=-[(d̃xk)t(Qxk+c)][(d̃xk)tQdxk].(iii)
The step sizeβLS in (3.55) is determined from the Armijo line search:(3.55)βLS=βk,
so that(3.56)Fμ(x̂k+1(βk))≤Fμ(x̂k)+βk∇Fμ(x̂k)Td̃x̂k,
where x̂k+1(βk)=(xk+1(βk),rk+1(βk),zk+1(βk))T=(xk+βkd̃xk,rk+βkd̃rk,zk+βkd̃zk)T; d̃x̂k=(d̃xk,d̃rk,d̃zk); so that d̃xk,d̃rk,d̃zkare defined in (3.29)–(3.31), respectively. In the proposed IPLS method, the Armijo search is also performed in the direction defined by primal variables. Starting at a point x̂k=(xk,rk,zk), the method searches for a new point x̂k+1, in direction d̃x̂k, so that the function Fμ defined in (3.3) decreases. This search calculates βLS, thereby ensuring the reduction of Fμ. To prevent oscillations in the iterative process, the initial choice for βk should not be too high, and to prevent the process from stopping prematurely, it should not be too low. Therefore, this value is generally adjusted to β0=1. If (3.44) is not satisfied, βk is updated using the following sequence:(3.57)βk+1=βkδ,k=0,1,…,n;withδ>1,usually set toδ=2.The following equation is used for the analysis of (3.56):
(3.58)∇Fμ(x̂k)=(∇Fμ(xk),∇Fμ(rk),∇Fμ(zk))T=(Qxk+c,-μkRk-1e,-μkZk-1e)T
with μk determined as in Section 3.6.The step size for dual variablesβDk is calculated by (3.59), without violating the nonnegativity requirements of the dual variables:(3.59)βDk=Min⁡{1,-αsikd̃sik,-αyikd̃yiksuch thatd̃sik;d̃yik<0}.The general principle used here to calculateαPk and βPk 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.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:(3.35)xk+1=xk+αPkdxk,(3.36)rk+1=rk+αPkdrk,(3.37)zk+1=zk+αPkdzk,(3.38)wk+1=wk+αDkdwk,(3.39)yk+1=yk+αDkdyk,(3.40)sk+1=sk+αDkdsk.The step sizeαPk for the primal variables is calculated through(3.41)αPk=Min⁡{αP,αQ,αLS},
where, for 0<α<1, the step sizes αP,αQ,αLS are determined as follows.(i)
The step size for primal variablesαP is obtained without violating the nonnegativity requirements of the primal variables:(3.42)αP=Min⁡{-αzikdzik,-αrikdriksuch thatdzik,drik<0}.(ii)
The step sizeαQ is the maximum step size possible without increasing the objective value:(3.43)αQ=-[(dxk)t(Qxk+c)][(dxk)tQdxk].(iii)
The step sizeαLS is determined as in (3.44) from the Fibonacci line search strategy, which is briefly summarized as follows:
(3.44)αLS=Min⁡αk⁡{Fμ(x̂k+1(αk))}=Min⁡αk⁡{Fμ(x̂k+αkdx̂k)},
where Fμ is defined in (3.3).For notation simplicity, in the summary of the Fibonacci search, the following identities are defined:x̂k+1(αLS)=(xk+1(αLS),rk+1(αLS),zk+1(αLS))T=(xk+αLSdxk,rk+αLSdrk,zk+αLSdzk)T, and dx̂k=(dxk,drk,dzk); so that dxkis defined in (3.23), drkis defined in (3.18), and dzkis defined in (3.17). Only primal variables are considered in the Fibonacci search algorithm used by the IPLS method. Starting at a point x̂k=(xk,rk,zk)T, the algorithm searches a new point x̂k+1 in direction dx̂k using the function Fμ defined in (3.3). The Fibonacci method calculates αLS so that the minimization of Fμ is ensured. αLS is determined considering the Fibonacci sequence. The initial value of αLS is set taking into account the interval of uncertainty [0,1].The step sizeαDk for the dual variables (3.38)–(3.40) is calculated through (3.45) for 0<α<1:(3.45)αDk=Min⁡{1,-αsikdsik,-αyikdyiksuch thatdsik;dyik<0}.


## 3.4.2. Corrector Step Size

In the corrector step, the calculation of step size is defined analogously, considering the variables from the following(3.46)xk+1=xk+βPkd̃xk,(3.47)rk+1=rk+βPkd̃rk,(3.48)zk+1=zk+βPkd̃zk,(3.49)wk+1=wk+βDkd̃wk,(3.50)yk+1=yk+βDkd̃yk,(3.51)sk+1=sk+βDkd̃yk,
where(3.52)βPk=Min⁡{βP,βQ,βLS},
while the step sizes βP,βQ,βLS are determined as follows.(i)
The step size for primal variablesβP is obtained without violating the nonnegativity requirements of the primal variables:(3.53)βP=Min⁡{1,-αzikd̃zik,-αrikd̃riksuch thatd̃zik,d̃rik<0}.(ii)
The step sizeβQ is the maximum step size possible without increasing the objective value:(3.54)βQ=-[(d̃xk)t(Qxk+c)][(d̃xk)tQdxk].(iii)
The step sizeβLS in (3.55) is determined from the Armijo line search:(3.55)βLS=βk,
so that(3.56)Fμ(x̂k+1(βk))≤Fμ(x̂k)+βk∇Fμ(x̂k)Td̃x̂k,
where x̂k+1(βk)=(xk+1(βk),rk+1(βk),zk+1(βk))T=(xk+βkd̃xk,rk+βkd̃rk,zk+βkd̃zk)T; d̃x̂k=(d̃xk,d̃rk,d̃zk); so that d̃xk,d̃rk,d̃zkare defined in (3.29)–(3.31), respectively. In the proposed IPLS method, the Armijo search is also performed in the direction defined by primal variables. Starting at a point x̂k=(xk,rk,zk), the method searches for a new point x̂k+1, in direction d̃x̂k, so that the function Fμ defined in (3.3) decreases. This search calculates βLS, thereby ensuring the reduction of Fμ. To prevent oscillations in the iterative process, the initial choice for βk should not be too high, and to prevent the process from stopping prematurely, it should not be too low. Therefore, this value is generally adjusted to β0=1. If (3.44) is not satisfied, βk is updated using the following sequence:(3.57)βk+1=βkδ,k=0,1,…,n;withδ>1,usually set toδ=2.The following equation is used for the analysis of (3.56):
(3.58)∇Fμ(x̂k)=(∇Fμ(xk),∇Fμ(rk),∇Fμ(zk))T=(Qxk+c,-μkRk-1e,-μkZk-1e)T
with μk determined as in Section 3.6.The step size for dual variablesβDk is calculated by (3.59), without violating the nonnegativity requirements of the dual variables:(3.59)βDk=Min⁡{1,-αsikd̃sik,-αyikd̃yiksuch thatd̃sik;d̃yik<0}.The general principle used here to calculateαPk and βPk 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 residualstk,uk, and also for the dual metric μk. Nevertheless, it is possible to use relative values for the metrics tk,uk, and μk 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:(3.60)‖tk‖‖b‖+1≤ε1(predictor step),(3.61)‖t̃k‖‖b‖+1≤ε1(corrector step);(ii)
dual feasibility:(3.62)‖uk‖‖Qxk+c‖+1≤ε2(predictor step),(3.63)‖ũk‖‖Qxk+c‖+1≤ε2(corrector step);(iii)
complementary slackness:(3.64)‖vk‖≤ε3,‖qk‖≤ε3(predictor step)(3.65)‖ṽk‖≤ε3,‖q̃k‖≤ε3(corrector step),
where ε1, ε2, and ε3 are sufficiently small positive numbers. Eventually, other criteria can be adopted according to the specific characteristics of each problem, as can be seen in [24, 32].


## 3.6. Update of the Barrier Parameter

According to [32], the barrier parameter is updated using an inner product that involves the primal variables rk and zk, and the dual variables sk and yk, respectively, as shown in the following equation(3.66)μk1=(rk)Tskn;μk2=(zk)Tykn,
so that μk is calculated using the following equation(3.67)μk=Min⁡{σμk1,σμk2}for a constant0<σ<1,
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).
Adjustk=0. Choose an arbitrary point:(x0;w0;z0,y0,r0,y0,s0)∈Ω0 and choose ε1, ε2, and ε3 as sufficiently small positive numbers.Step 2 (intermediate calculations—predictor).
Calculate:tk;gk;fk;ok;qk;vk using (3.10), μk using (3.67) and matrix θ using (3.21).Step 3 (finding directions of translation—predictor).
Calculate search directionsdxk, dwk, dzk, drk, dyk, and dsk for the predictor step using (3.17)–(3.24).Step 4 (computing step size—predictor).
Calculate Fibonacci step sizesαPk and αDk using (3.41)–(3.44).Step 5 (moving to a new solution—predictor).
Updatexk+1; wk+1; sk+1; yk+1; zk+1;rk+1 obtained from the predictor step, according to (3.35)–(3.44).Step 6 (checking for optimality—predictor).
If the criteria defined in (3.60), (3.62), and (3.64) are satisfied, stop. The solution is optimal. Otherwise, go on to the following step.Step 7 (intermediate calculations—corrector).
Calculate:t̃k; g̃k; f̃k; õk; q̃k; ṽk using (3.27), μk using (3.67) and matrix θ using (3.21).Step 8 (finding directions of translation—corrector).
Calculated̃xk, d̃wk, d̃zk, d̃rk, d̃yk, and d̃sk, for the corrector step using (3.28)–(3.33).Step 9 (computing step size—corrector).
Calculate the Armijo step sizeβPk and βDk using (3.52)–(3.58).Step 10 (moving to a new solution—corrector).
Updatexk+1; wk+1; sk+1; yk+1; zk+1; rk+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.
Adjustk=k+1 and return to Step 2.The predictor Steps2 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. Tables1, 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.
UnitPimin (MW)Pimax  (MW)ai ($/MW2)bi ($/MW)ci ($)11006000,0015627,925612502000,0048207,977831004000,0019407,85310Tables2 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.
Power outputHGACOEGAPCPDIPLSP1 (MW)470,8421344,7295393,1698393,1698P2 (MW)109,4012193,9445122,2264122,2264P3 (MW)269,7567311,3260334,6038334,6038ΣPi  (MW)850,0000850,0000850,000850,000Table 3
Characteristics of the system with 3 generators.
ResultsGeneration costs ($)HGA8.212,73COEGA8.223,86PCPD8.194,36IPLS8.194,36The 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 Table10. 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 Table1. Parameters ai,bi,ci stand for the coefficients of the cost functions for the generators, while Pimin,Pimax represent minimum and maximum power output capabilities, respectively, for generating unit i.The following values are adopted to initialise the method:(4.1)x0=(450,100,300);w0=(0,0,0);y0=(0,0,0).The parameters related to the system’s total demand are set atD=850MW and the active power losses are neglected. The values adopted for ε1,ε2, and ε3 are 10-8.Table2 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

Table4 describes the main characteristics of the system containing 6 generators.Table 4
Characteristics of the system with 6 generators.
UnitPiMin (MW)PiMax  (MW)ai ($/MW2)bi ($/MW)ci ($)1101250,1524738,53973756,798862101500,1058746,15916451,32513352250,0280340,396551049,99774352100,0354638,305331243,531151303250,0211136,327821658,569661253150,0179938,240411356,6592Table 5
Comparison of the power generation outputs obtained for the system with 6 generators.
Generated powerHGAPCPDIPLSP1 (MW)20,136717,3659617,36596P2 (MW)14,864510,0000010,00000P3 (MW)72,400861,34066761,340667P4 (MW)72,400877,9748777,97487P5 (MW)180,0617177,81828177,81828P6 (MW)139,5865155,500216155,500216ΣPi (MW)500,0000500,0000500,0000Table 6
Values of the objective function for the system with 6 generators.
ResultsGeneration cost ($)HGA27.037,29PCPD27.003,50IPLS27.003,50As in the previous case study, the parametersai,bi,ci stand for the coefficients of the cost functions for the generating unit I, while Pimin,Pimax represent minimum and maximum power output capabilities, respectively, for the generating unit i.The following values are adopted to initialise the method:(4.2)x0=(20,30,75,75,145,155);w0=(0,0,0,0,0,0);y0=(1,1,1,1,1,1).The parameters related to the system’s total demand are set atD=500MW and the active power losses are neglected. The values adopted for ε1,ε2, and ε3 are 10-8.Table5 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

Table7 describes the main characteristics of the system containing 13 generators.Table 7
Characteristics of the system with 13 generators.
UnitPiMin (MW)PiMax (MW)ai ($/MW2)bi  ($/MW)ci ($)106800,000288,1550203600,000568,1309303600,000568,13074601800,003247,742405601800,003247,742406601800,003247,742407601800,003247,742408601800,003247,742409601800,003247,7424010401200,002848,612611401200,002848,612612551200,002848,612613551200,002848,6126Table 8
Comparison of the power generation outputs obtained for the system with 13 generators.
Power outputHGACOEGACAPCPDIPLSP1651,145735,626679,2551680,0000680,0000P2319,982337,495359,8672360,0000360,0000P3320,463292,625357,2368360,0000360,0000P4137,776146,713154,8137155,0000155,0000P5156,688177,346158,0946155,0000155,0000P6147,007131,552155,8520155,0000155,0000P7159,165154,197146,1697155,0000155,0000P8145,378159,550146,8364155,0000155,0000P9151,551167,339168,7979155,0000155,0000P1082,25960,67740,018140,000040,0000P1186,32074,68140,000040,000040,0000P1282,89356,53755,017555,000055,0000P1379,36825,65555,048855,000055,0000ΣPi  (MW)2.520,002.520,002.520,002.520,002.520,00Table 9
Values of the objective function for the system with 13 generators.
ResultsGeneration cost ($)HGA24.111,69COEGA24.072,03HAGA24.052,34CA24.052,10PCPD24.050,08IPLS24.050,08Table 10
Number of iterations for the PCPD and IPLS methods.
Simulated systemPCPDIPLS3 generating units26126 generating units221513 generating units2317The parameters in this table are analogous to those described in the preceding case studies. The following values are adopted to initialise the method:(4.3)x0=(660,320,330,150,140,155,165,150,140,70,80,75,85);w0=(0,0,0,0,0,0,0,0,0,0,0,0,0);y0=(30,30,30,30,30,30,30,30,30,30,30,30,30).The parameters related to the system’s total demand are set atD=2520MW and the active power losses are neglected. The values adopted for ε1,ε2, and ε3 are 10-8.Table8 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, Table10 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 Table10.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:(4.4)MinimizeCTsubject to∑i=1nPi=DPi∈S,
where S={Pi∈ℜ⊥Pi-Pimin≥0andPimax-Pi≥0,i=1…,n}, and CT = ∑i=1n(Ci(Pi) = aiPi2+biPi+ci).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:(4.5)L(Pi,w,si,yi)=CT(Pi)+w(D-∑i=1nPi)+∑i=1nsi(Pimin⁡-Pi)+∑i=1nyi(Pi-Pimax⁡),
where w∈ℜ,  si∈ℜ, and yi∈ℜ,  i=1,…,n, are the associated Lagrange multipliers. The KKT conditions associated with the Lagrangian function (4.5) are expressed through(4.6)D-∑i=1nPi=0,(4.7)Pi-Pimin⁡≥0;si(Pi-Pimin⁡)=0;i=1,…,n,(4.8)Pimax⁡-Pi≥0;yi(Pimax⁡-Pi)=0;i=1,…,n,(4.9)λi(Pi)-w-si+yi=0;i=1,…,n,
where(4.10)λi(Pi)=∂Ci(Pi)∂Pi=2aiPi+bi;i=1,…,n.In the literature on power systems, the Lagrange multiplierw is called energy price (i.e., the price of 1 MWh), while λi(Pi) is called incremental or marginal cost of generating unit i. Using expression (4.9), the energy price can be calculated by(4.11)w=λi(Pi)-si+yi;i=1,…,n.


### 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 thatsi*=0, yi*=0; i=1,…,n. Thus, from (4.11), one reaches(4.12)w=λi*(Pi*);i=1,…,n,
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, si*=0 for i=1,…,n, yj*>0 for j∈Ωmax, and yj*=0for j∉Ωmax. Starting from (4.11), one reaches(4.13)w=λj*(Pj*)+yj*;j∈Ωmax⁡,w=λj*(Pj*);j∉Ωmax⁡.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 yj*>0;j∈Ω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, yi*=0 for i=1,…,n, sj*>0 for j∈Ωmin, and sj*=0 for j∉Ωmin. Thus, using (4.11), one reaches(4.14)w=λj*(Pj*)-sj*;j∈Ωmin⁡,w=λj*(Pj*);j∉Ωmin⁡.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 sj*>0;j∈Ω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

Table2 describes the optimal dispatch, Pi*,i=1,…,n, calculated by the IPLS method for the system with 3 generating units. As this table indicates, no generating unit has reached its 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 si* and yi*;i=1,…,n are all zero. Also, the energy price is equivalent to the marginal costs, which are all equal, as expected. These results are in accordance with the theory described in Section 4.7.1.Table 11
Incremental analysis for the 3-generator system.
Unitλi*si*yi*w*19,14826009,1482629,148260039,1482600


#### 4.8.2. Case 2: Power System with 6 Generating Units

Table5 shows the optimal dispatch, Pi*,i=1,…,n, calculated by the IPLS method for the system with 6 generating units. As can be seen in this table, the generating unit 2 has reached its lower limit, so the optimal solution has one active constraint. This situation coincides with the one described in Section 4.7.3. Table 12 lists the energy price, marginal costs, and Lagrange multipliers calculated by the IPLS method for this case. As described in Section 4.7.3, the values obtained by the IPLS method for si* and yi*;i=1,…,n are all zero except for the value of s2*, 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.
Unitλi*si*yi*w143,84490043,8449248,27654,431640343,844900443,844900543,844900643,844900


#### 4.8.3. Case 3: Power System with 13 Generating Units

Table8 shows the optimal dispatch, Pi*,i=1,…,n, calculated by the IPLS method for the system with 13 generating units. Table 13 lists the energy price, marginal costs, and Lagrange multipliers calculated by the IPLS method. As the latter table indicates, generated units 1, 2, and 3 have reached their upper limits, and generated units 10, 11, 12, and 13 have reached their lower limits, indicating that there are 6 active constraints in the optimal solution calculated by the IPLS method. This situation coincides with the one described in Section 4.7.4. As described previously, the values obtained by the IPLS method for si* andyi*; i=1,…,n are zero except for the values of y1*,y2*,y3*, which are associated with upper output capabilities, and also for s10*,s11*,s12*,s13*, 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.
Unitλi*si*yi*w18,480700,26358,744328,503100,241138,503100,241148,74430058,74430068,74430078,74430088,74430098,744300108,82720,08280118,82720,08280128,91240,16800138,91240,16800As expected, in all the cases analysed in Tables11, 12, and 13, the value of the generation cost w=λi*-si*+yi* was determined univocally, which proves the equal incremental cost criterion [29].The analysed results indicate that the dispatchesPi*,i=1,…,n calculated, respectively, in Tables 2, 5, and 8 and the results for λi*, w, si*, and yi* determined, respectively, in Tables 11, 12, and 13, satisfy the complementary slackness conditions set forth in (4.7) and (4.8).


## 4.1. Power System with 3 Generators

The main characteristics of the system containing 3 generators are described in Table1. Parameters ai,bi,ci stand for the coefficients of the cost functions for the generators, while Pimin,Pimax represent minimum and maximum power output capabilities, respectively, for generating unit i.The following values are adopted to initialise the method:(4.1)x0=(450,100,300);w0=(0,0,0);y0=(0,0,0).The parameters related to the system’s total demand are set atD=850MW and the active power losses are neglected. The values adopted for ε1,ε2, and ε3 are 10-8.Table2 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

Table4 describes the main characteristics of the system containing 6 generators.Table 4
Characteristics of the system with 6 generators.
UnitPiMin (MW)PiMax  (MW)ai ($/MW2)bi ($/MW)ci ($)1101250,1524738,53973756,798862101500,1058746,15916451,32513352250,0280340,396551049,99774352100,0354638,305331243,531151303250,0211136,327821658,569661253150,0179938,240411356,6592Table 5
Comparison of the power generation outputs obtained for the system with 6 generators.
Generated powerHGAPCPDIPLSP1 (MW)20,136717,3659617,36596P2 (MW)14,864510,0000010,00000P3 (MW)72,400861,34066761,340667P4 (MW)72,400877,9748777,97487P5 (MW)180,0617177,81828177,81828P6 (MW)139,5865155,500216155,500216ΣPi (MW)500,0000500,0000500,0000Table 6
Values of the objective function for the system with 6 generators.
ResultsGeneration cost ($)HGA27.037,29PCPD27.003,50IPLS27.003,50As in the previous case study, the parametersai,bi,ci stand for the coefficients of the cost functions for the generating unit I, while Pimin,Pimax represent minimum and maximum power output capabilities, respectively, for the generating unit i.The following values are adopted to initialise the method:(4.2)x0=(20,30,75,75,145,155);w0=(0,0,0,0,0,0);y0=(1,1,1,1,1,1).The parameters related to the system’s total demand are set atD=500MW and the active power losses are neglected. The values adopted for ε1,ε2, and ε3 are 10-8.Table5 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

Table7 describes the main characteristics of the system containing 13 generators.Table 7
Characteristics of the system with 13 generators.
UnitPiMin (MW)PiMax (MW)ai ($/MW2)bi  ($/MW)ci ($)106800,000288,1550203600,000568,1309303600,000568,13074601800,003247,742405601800,003247,742406601800,003247,742407601800,003247,742408601800,003247,742409601800,003247,7424010401200,002848,612611401200,002848,612612551200,002848,612613551200,002848,6126Table 8
Comparison of the power generation outputs obtained for the system with 13 generators.
Power outputHGACOEGACAPCPDIPLSP1651,145735,626679,2551680,0000680,0000P2319,982337,495359,8672360,0000360,0000P3320,463292,625357,2368360,0000360,0000P4137,776146,713154,8137155,0000155,0000P5156,688177,346158,0946155,0000155,0000P6147,007131,552155,8520155,0000155,0000P7159,165154,197146,1697155,0000155,0000P8145,378159,550146,8364155,0000155,0000P9151,551167,339168,7979155,0000155,0000P1082,25960,67740,018140,000040,0000P1186,32074,68140,000040,000040,0000P1282,89356,53755,017555,000055,0000P1379,36825,65555,048855,000055,0000ΣPi  (MW)2.520,002.520,002.520,002.520,002.520,00Table 9
Values of the objective function for the system with 13 generators.
ResultsGeneration cost ($)HGA24.111,69COEGA24.072,03HAGA24.052,34CA24.052,10PCPD24.050,08IPLS24.050,08Table 10
Number of iterations for the PCPD and IPLS methods.
Simulated systemPCPDIPLS3 generating units26126 generating units221513 generating units2317The parameters in this table are analogous to those described in the preceding case studies. The following values are adopted to initialise the method:(4.3)x0=(660,320,330,150,140,155,165,150,140,70,80,75,85);w0=(0,0,0,0,0,0,0,0,0,0,0,0,0);y0=(30,30,30,30,30,30,30,30,30,30,30,30,30).The parameters related to the system’s total demand are set atD=2520MW and the active power losses are neglected. The values adopted for ε1,ε2, and ε3 are 10-8.Table8 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, Table10 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 Table10.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:(4.4)MinimizeCTsubject to∑i=1nPi=DPi∈S,
where S={Pi∈ℜ⊥Pi-Pimin≥0andPimax-Pi≥0,i=1…,n}, and CT = ∑i=1n(Ci(Pi) = aiPi2+biPi+ci).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:(4.5)L(Pi,w,si,yi)=CT(Pi)+w(D-∑i=1nPi)+∑i=1nsi(Pimin⁡-Pi)+∑i=1nyi(Pi-Pimax⁡),
where w∈ℜ,  si∈ℜ, and yi∈ℜ,  i=1,…,n, are the associated Lagrange multipliers. The KKT conditions associated with the Lagrangian function (4.5) are expressed through(4.6)D-∑i=1nPi=0,(4.7)Pi-Pimin⁡≥0;si(Pi-Pimin⁡)=0;i=1,…,n,(4.8)Pimax⁡-Pi≥0;yi(Pimax⁡-Pi)=0;i=1,…,n,(4.9)λi(Pi)-w-si+yi=0;i=1,…,n,
where(4.10)λi(Pi)=∂Ci(Pi)∂Pi=2aiPi+bi;i=1,…,n.In the literature on power systems, the Lagrange multiplierw is called energy price (i.e., the price of 1 MWh), while λi(Pi) is called incremental or marginal cost of generating unit i. Using expression (4.9), the energy price can be calculated by(4.11)w=λi(Pi)-si+yi;i=1,…,n.


## 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 thatsi*=0, yi*=0; i=1,…,n. Thus, from (4.11), one reaches(4.12)w=λi*(Pi*);i=1,…,n,
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, si*=0 for i=1,…,n, yj*>0 for j∈Ωmax, and yj*=0for j∉Ωmax. Starting from (4.11), one reaches(4.13)w=λj*(Pj*)+yj*;j∈Ωmax⁡,w=λj*(Pj*);j∉Ωmax⁡.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 yj*>0;j∈Ω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, yi*=0 for i=1,…,n, sj*>0 for j∈Ωmin, and sj*=0 for j∉Ωmin. Thus, using (4.11), one reaches(4.14)w=λj*(Pj*)-sj*;j∈Ωmin⁡,w=λj*(Pj*);j∉Ωmin⁡.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 sj*>0;j∈Ω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.7.1. No Active Constraint

In this case there is no active constraint in the optimal solution of problem (4.4), so thatsi*=0, yi*=0; i=1,…,n. Thus, from (4.11), one reaches(4.12)w=λi*(Pi*);i=1,…,n,
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, si*=0 for i=1,…,n, yj*>0 for j∈Ωmax, and yj*=0for j∉Ωmax. Starting from (4.11), one reaches(4.13)w=λj*(Pj*)+yj*;j∈Ωmax⁡,w=λj*(Pj*);j∉Ωmax⁡.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 yj*>0;j∈Ω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, yi*=0 for i=1,…,n, sj*>0 for j∈Ωmin, and sj*=0 for j∉Ωmin. Thus, using (4.11), one reaches(4.14)w=λj*(Pj*)-sj*;j∈Ωmin⁡,w=λj*(Pj*);j∉Ωmin⁡.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 sj*>0;j∈Ω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

Table2 describes the optimal dispatch, Pi*,i=1,…,n, calculated by the IPLS method for the system with 3 generating units. As this table indicates, no generating unit has reached its 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 si* and yi*;i=1,…,n are all zero. Also, the energy price is equivalent to the marginal costs, which are all equal, as expected. These results are in accordance with the theory described in Section 4.7.1.Table 11
Incremental analysis for the 3-generator system.
Unitλi*si*yi*w*19,14826009,1482629,148260039,1482600


### 4.8.2. Case 2: Power System with 6 Generating Units

Table5 shows the optimal dispatch, Pi*,i=1,…,n, calculated by the IPLS method for the system with 6 generating units. As can be seen in this table, the generating unit 2 has reached its lower limit, so the optimal solution has one active constraint. This situation coincides with the one described in Section 4.7.3. Table 12 lists the energy price, marginal costs, and Lagrange multipliers calculated by the IPLS method for this case. As described in Section 4.7.3, the values obtained by the IPLS method for si* and yi*;i=1,…,n are all zero except for the value of s2*, 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.
Unitλi*si*yi*w143,84490043,8449248,27654,431640343,844900443,844900543,844900643,844900


### 4.8.3. Case 3: Power System with 13 Generating Units

Table8 shows the optimal dispatch, Pi*,i=1,…,n, calculated by the IPLS method for the system with 13 generating units. Table 13 lists the energy price, marginal costs, and Lagrange multipliers calculated by the IPLS method. As the latter table indicates, generated units 1, 2, and 3 have reached their upper limits, and generated units 10, 11, 12, and 13 have reached their lower limits, indicating that there are 6 active constraints in the optimal solution calculated by the IPLS method. This situation coincides with the one described in Section 4.7.4. As described previously, the values obtained by the IPLS method for si* andyi*; i=1,…,n are zero except for the values of y1*,y2*,y3*, which are associated with upper output capabilities, and also for s10*,s11*,s12*,s13*, 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.
Unitλi*si*yi*w18,480700,26358,744328,503100,241138,503100,241148,74430058,74430068,74430078,74430088,74430098,744300108,82720,08280118,82720,08280128,91240,16800138,91240,16800As expected, in all the cases analysed in Tables11, 12, and 13, the value of the generation cost w=λi*-si*+yi* was determined univocally, which proves the equal incremental cost criterion [29].The analysed results indicate that the dispatchesPi*,i=1,…,n calculated, respectively, in Tables 2, 5, and 8 and the results for λi*, w, si*, and yi* determined, respectively, in Tables 11, 12, and 13, satisfy the complementary slackness conditions set forth in (4.7) and (4.8).


## 4.8.1. Case 1: Power System with 3 Generating Units

Table2 describes the optimal dispatch, Pi*,i=1,…,n, calculated by the IPLS method for the system with 3 generating units. As this table indicates, no generating unit has reached its 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 si* and yi*;i=1,…,n are all zero. Also, the energy price is equivalent to the marginal costs, which are all equal, as expected. These results are in accordance with the theory described in Section 4.7.1.Table 11
Incremental analysis for the 3-generator system.
Unitλi*si*yi*w*19,14826009,1482629,148260039,1482600


## 4.8.2. Case 2: Power System with 6 Generating Units

Table5 shows the optimal dispatch, Pi*,i=1,…,n, calculated by the IPLS method for the system with 6 generating units. As can be seen in this table, the generating unit 2 has reached its lower limit, so the optimal solution has one active constraint. This situation coincides with the one described in Section 4.7.3. Table 12 lists the energy price, marginal costs, and Lagrange multipliers calculated by the IPLS method for this case. As described in Section 4.7.3, the values obtained by the IPLS method for si* and yi*;i=1,…,n are all zero except for the value of s2*, 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.
Unitλi*si*yi*w143,84490043,8449248,27654,431640343,844900443,844900543,844900643,844900


## 4.8.3. Case 3: Power System with 13 Generating Units

Table8 shows the optimal dispatch, Pi*,i=1,…,n, calculated by the IPLS method for the system with 13 generating units. Table 13 lists the energy price, marginal costs, and Lagrange multipliers calculated by the IPLS method. As the latter table indicates, generated units 1, 2, and 3 have reached their upper limits, and generated units 10, 11, 12, and 13 have reached their lower limits, indicating that there are 6 active constraints in the optimal solution calculated by the IPLS method. This situation coincides with the one described in Section 4.7.4. As described previously, the values obtained by the IPLS method for si* andyi*; i=1,…,n are zero except for the values of y1*,y2*,y3*, which are associated with upper output capabilities, and also for s10*,s11*,s12*,s13*, 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.
Unitλi*si*yi*w18,480700,26358,744328,503100,241138,503100,241148,74430058,74430068,74430078,74430088,74430098,744300108,82720,08280118,82720,08280128,91240,16800138,91240,16800As expected, in all the cases analysed in Tables11, 12, and 13, the value of the generation cost w=λi*-si*+yi* was determined univocally, which proves the equal incremental cost criterion [29].The analysed results indicate that the dispatchesPi*,i=1,…,n calculated, respectively, in Tables 2, 5, and 8 and the results for λi*, w, si*, and yi* 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.



---

*Source: 376546-2012-06-04.xml*