- Email: [email protected]
Modified Newton method for reactive dispatching
Electrical Power and Energy Systems 24 (2002) 815±819
www.elsevier.com/locate/ijepes
Modi®ed Newton method for reactive dispatching G.R.M. da Costa Department of Electrical Engineering, Escola de Engenharia de SaÄo Carlos, USP, 13560-970 SaÄo Carlos SP, Brazil Received 22 December 2000; received in revised form 11 October 2001; accepted 13 December 2001
Abstract This paper describes a new approach to the optimal reactive dispatch problem, based on an augmented Lagrangian function of the original problem. The Karush±Kuhn±Tucker (KKT) optimality conditions are solved by the modi®ed Newton method. The second-order information in the original system of equations is approximated and the ®rst-order information is kept intact The proposed method requires less computer memory than those algorithms currently available. The effectiveness of the proposed approach has been examined by solving the IEEE 30bus and the BRAZILIAN 810-bus systems. q 2002 Elsevier Science Ltd. All rights reserved. Keywords: Power systems; Optimal power ¯ow; Optimization; Multiplier method; Augmented Lagrangian
1. Introduction The optimal reactive dispatch (ORD) is a particular case of the optimal power ¯ow (OPF) in which active controls are ®xed and control optimization only relates to reactive power, e.g. voltage magnitude for generators, synchronous condensers and static VAr, taps and switchable banks. The OPF problem was proposed by Carpentier in the early 1960s based on the economic dispatch problem [1]. Since then, many papers have been written in the attempt to solve it. Gradient techniques were the ®rst approach applied to the problem [2]. Slow convergence limited their use. Sasson et al. [3] proposed a method in which the objective function is penalized with both equality and inequality constraints. The Hessian matrix was used for the ®rst time. Sun et al. [4] introduced into OPF research a combination of the Newton's method with Lagrange multipliers and penalty functions. The major dif®culty turned out to be the identi®cation of binding constraints. The Dual-Newton approach [5] handles equality constraints with Lagrange multipliers and inequality constraints with both penalty factors and Lagrange multipliers. The binding constraints need not be identi®ed. In Ref. [6] a primal±dual (PD) logarithmic barrier algorithm is applied directly to a nonlinear OPF problem and the Karush±Kuhn±Tucker (KKT) conditions are solved by Newton's method. A crucial step in the algorithm is the choice of the barrier parameter. Costa et al. [7] showed a new approach that improves the performance of the Newton's method proposed by Sun et al. [4]. Up to the E-mail address: [email protected] (G.R.M. da Costa).
present there is no single robust, reliable and fast approach that ful®lls the requirements of the EMS. Huneault and Galiana [8] have made an extensive study of the approaches used to solve the OPF problem. It is evident that all these solutions have followed close behind advances in numerical optimization, and where possible were computed by techniques already well established, which use the incidence matrix structure. This structure has the advantage that computer packages can be applied to it. In this paper, we propose a new approach to the solution of the ORD problem. An augmented Lagrangian is associated with the original problem and the KKT optimality conditions are solved by the modi®ed Newton method. The constraints include the reactive power limits of generators, the limits on the voltage at the load buses, and the limits of the control variables, i.e. the transformer tap positions, generator terminal voltages, switchable reactive power sources and limits on transmission line ¯ow. All variables are treated as continuous. The paper is organized as follows: ®rst, the problem of ORD is shown. Next, the modi®ed Newton and the augmented Lagrangian approaches are discussed. The implementation of this new approach is described. The results obtained testing on two systems (30 and 810 buses) are reported. Finally, some conclusions are drawn.
2. Description of the ORD problem Minimization of real power losses in a system forms the basis of the reactive power optimization problem. The ORD
0142-0615/02/$ - see front matter q 2002 Elsevier Science Ltd. All rights reserved. PII: S 0142-061 5(02)00013-3
816
G.R.M. da Costa / Electrical Power and Energy Systems 24 (2002) 815±819
problem can be presented as:
constraints of the MORD model is de®ned by
minimize f x
L f x 1
i1
such that gi x 0 i 1; 2; ¼; m , n hj x # 0 j 1; 2; ¼; p
m X
(1)
where x [ R n is the vector of state variable, f(x) the real power losses of the transmission, g x 0 the power ¯ow equations, and h x # 0 is the limits on state variable and functional constraints of the power system. The state variable vector x represents the voltage magnitude, angles and LTCs taps. The objective function, f(x), represents the active power losses in the transmission. This function is nonseparable and permits no simpli®cations. The equality constraints, g x 0; represent the power ¯ow equations for scheduled load and generation. The inequality constraints, h x # 0; represent the state variables such as voltage magnitude limits and functional constraints of the power ¯ow, e.g. limit of active and reactive power ¯ows in the transmission lines and transformers, limits of reactive power injections for reactive control buses and active power injection for the slack bus. This is a typical nonlinear and nonconvex problem.
1
li gi x
p X v mj hj x 1 ej 1 hj x 1 ej 2 2 j1
where l i, m j are the Lagrange multipliers (related to the equality and inequality constraints, respectively). The minimum of L with respect to the slack variables (ej) is calculated from the following stationary conditions: 2L 0 2ej
4
that is:
mj 1 v h x 1 ej 0 j 1; ¼; p These conditions are satis®ed by mj ej max 0; 2 2 hj x v
p v X h x 1 ej 2 2 j1
such that gi x 0
6
MORD : p v X r 2 x 2 j1 j
7
such that gi x 0
When the problem (1) is analyzed in terms of its ordinary Lagrangian function, there is usually a duality gap, unless the objective and constraint functions are convex. The gap can be removed by passing to an augmented Lagrangian function that involves quadratic penalty-like terms. The modi®ed dual problem then consists of maximizing a concave function of the Lagrange multipliers and an additional variable, which is a penalty parameter. The problem becomes a sequential unconstrained minimization problem solved by Newton-based methods. The following modi®ed optimal reactive dispatch (MORD) model, related with the ORD, is obtained then as: MORD : minimize f x 1
5
Finally, the expression above is substituted in the MORD model, giving the following:
minimize f x 1 3. The modi®ed Newton method applied to an augmented Lagrangian function
3
2
hj x 1 ej 0 ej $ 0 where ej is the slack variable related to inequality constraint, hj(x) and v is the penalty factor. The classical Lagrangian corresponding to the equality
rj x # 0 where
8 > > < hj x;
mj v rj x m m > j > : 2 ; if hj x # 2 j v v if hj x $ 2
8
The Lagrangian related to the MORD model, as formulated in Eq. (7), is given by: p m X X v li gi x 1 rj x mj 1 rj x 9 L f x 1 2 i1 j1 The Augmented Lagrangian relating to the MORD model is de®ned by Eq. (9). The solution of the MORD model will be sought by an operational strategy de®ned on conditions imposed on the function L. The process involves the determination of values for x and l , such that the conditions of optimality are satis®ed. Thus: 2L 0; 2x
2L 0 2l
10
This system of equations was solved using the modi®ed Newton method in which Dx and Dl were calculated.
G.R.M. da Costa / Electrical Power and Energy Systems 24 (2002) 815±819
Thus, the Eq. (10) can be represented as: "
B
JT
J
0
#"
Dx Dl
#
" 2
7xL
slowly, to in®nity:
# 11
7 lL
where B is an approximation to the Hessian matrix, i.e. B is a diagonal matrix of second partial derivatives of L with respect to x. Each element of B is a linear combination of second partial derivatives of power ¯ow equations, its constraints and objective function. J is a matrix of ®rst partial derivatives of power ¯ow equations with respect to x. The second-order information in the original system of equations is approximated and the ®rst-order information is kept intact. There are various ways to choose the approximation B [10]. Here we assume B equal to the diagonal of the Hessian matrix, updated to each iteration. One advantage of this method is that B can be taken to be positive de®nite even when the Hessian matrix, is not as given in Ref. [10]. For conciseness, Eq. (11) can be written as: WB Dz 2b
12
where WB is the Lagrangian matrix, Dz the vector of corrections formed by Dx and Dl , and b is the gradient of L with respect to x and l . 3.1. Updating the multipliers and penalty factors The gradient of the augmented Lagrangian at the minimum point x p is: 2L 2f xp 2g xp 1l 2x 2x 2x 8 2h xp > > < m 1 vh xp 2x 1 > > :0
m v m p if h x # 2 v
if h xp $ 2
mk11 j
> > > :0
if hj x $ 2
mkj v
mkj if hj x # 2 v
vk11 bvk
15
where b is the penalty increasing factor (b . 1). The ill-conditioning usually associated with the penalty function approach can be controlled. The method treats the equality and inequality constraints in a uni®ed manner, and avoids a search for binding constraints. 3.2. Algorithm The algorithm proposed to solve the ORD problem is an iterative process consisting of the following steps: (i) Make starting estimates for z x; l and m x can be the same as the initial values for a power ¯ow l m 0 or any reasonable guess (ii) Evaluate b as a function of z (iii) If KKT conditions are satis®ed the problem is solved, otherwise (iv) Update m (Eq.(14)) and v (Eq.(15)) (v) Evaluate the matrix WB as a function of z (vi) Solve the system WB Dz 2b for Dz (vii) Update z by Dz (viii) Return to step (ii). The initial values allocated to the variables x can be in the infeasible region of the problem or can be a power ¯ow solution. In the evaluation of matrix WB in step (v), the updated values of m and v must be used. For suf®ciently large v, WB is positive de®nite, thus a problem can be convexi®ed (at least locally).
4. Implementation 13
The terms m 1 vh x can be seen as the Lagrange multipliers that satisfy the optimality condition of the problem (1), [9]. Therefore, a good approximation to the unknown m is: 8 > k k11 > > < mj 1 vhj x ;
817
Most of the work in the algorithm is in the solution of system (11). The Lagrangian matrix, WB, that results from the linear approximation of the KKT conditions, has a structure that facilitates the application of sparsity of techniques. The three-bus system of Fig. 1 is used as an example to show the structure of WB (Eq. (16)). Each symbol indicates a nonzero element and its coordinates indicate a speci®c derivative. For example, B with coordinates u2 ; u2 is 22 L=2u22 and J with coordinates V2 ; lP3 is 22 L=2V2 2lP3 2DP3 =2V2 : Only the horizontal
14
The penalty parameter v may also be adjusted during the process. As in ordinary penalty function methods, the sequences of v increase towards a ®nite value, or tend,
Fig. 1. Three-bus electrical system.
818
G.R.M. da Costa / Electrical Power and Energy Systems 24 (2002) 815±819
Table 1 Speci®cations of tests cases
Table 3 Optimization summary for 30-bus system
Systems
Buses
Lines
Q-controls
LTC-controls
30-Bus 810-Bus
30 810
41 1340
5 131
4 205
coordinates are shown: v1 2 B 6 6 6 6 6 6 6 6 WB 6 6 6 6 6 6 6 6J 6 6 6 4
u2
J
u3
v2
v3
lP2
lP3
J B B B B J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
lQ2 3 J 7 J7 7 7 J7 7 7 7 J7 7 7 J7 7 7 7 7 7 7 7 5
Iteration
Active power loss (MW)
Mismatch uDPuMax (MW)
Mismatch uDQuMax (MVAr)
1 2 3 4 5 6 7 8 9 10 11
18.76 18.38 17.80 17.41 17.40 17.33 17.26 17.25 17.23 17.23 17.22
0.1763 0.0091 0.0220 0.0163 0.0083 0.0023 0.0001 0.0000 0.0000 0.0000 0.0000
0.1411 0.0956 0.0867 0.0225 0.0237 0.0381 0.0034 0.0024 0.0006 0.0002 0.0001
16
This matrix is sparse and symmetric. It needs to compute and store only half of LU factorization due to symmetry. The matrix structure is constant through iterations, the ordering and symbolic analysis are done only once to create a static data structure. Thus the numerical factorization is carried out ef®ciently at every iteration.
multipliers related to the equality and inequality constraints are, respectively, l 0 and m 0: Table 2 presents the sensitivity of the optimization process to the penalty factors, v, of the voltage magnitude. On the basis of these results, the initial penalties were de®ned as v 500 for voltage magnitude and ratio tap while the other inequality constraints were v 1:0: All penalty increasing factors were de®ned as b 1:2: The process converged in 11 iterations. CPU time for this run was 1.14 s. The amount of reactive power generation was 147.16 MVAr, with a total active power loss of 17.22 MW.
5. Test results Tests were done to verify the ef®ciency of the proposed approach. The algorithm was implemented in FORTRAN, using double precision arithmetic on a SUN SPARC 20 workstation, in the Power Systems Optimization Laboratory of EESC, USP. The cases studied were the minimization of active power losses in transmission in the IEEE 30-bus and the BRAZILIAN SOUTH/SOUTHEAST 810-bus systems. Table 1 lists the speci®cations of these test cases.
Fig. 2. Mismatches for the 810-bus system.
5.1. 30-Bus system This test was accomplished with the following initial conditions: Vk 1.0 pu (100 MVA base) and uk 0:0 for k 1; ¼; 30; and ti 1:0 for i 1; ¼; 4: The Lagrange Table 2 Tests for 30-bus system Tests
b
v
CPU time (s)
Iterations
1 2 3 4
1.2 1.1 1.2 1.1
100 300 500 1000
1.71 1.63 1.14 Diverges
19 16 11 Diverges
Fig. 3. Objective function for 810-bus system.
G.R.M. da Costa / Electrical Power and Energy Systems 24 (2002) 815±819 Table 4 Comparative tests on 810-bus system Algorithm
Losses (MW)
Reactive (MVAr)
Iterations
CPU Time (s)
ASP PDBL PD MN
1690.01 1702.17 1689.42 1693.03
2498.73 2510.11 2496.89 2499.70
19 39 28 46
217.27 270.58 259.95 204.31
The optimization process for this case is summarized in Table 3. 5.2. 810-Bus system This test was accomplished with the following initial conditions: x t; V; u unconverged network solution. The Lagrange multipliers for the equality and inequality constraints were, respectively, l 0 and m 0: The initial penalties were de®ned as v 2000 for voltage magnitude and ratio tap and the other inequality constraints were v 1:0: All penalty increasing factors were de®ned as b 1:04: The process converged in 46 iterations to a speci®ed tolerance j 10 24 : CPU time for this run was 204.31 s. The amount of reactive power generation was 2499.70 MVAr, with a total active power loss of 1693.03 MW. The load ¯ow does not keep the voltage magnitude within desirable limits on the load buses, with a total of 2839.69 MVAr of reactive power generation and 1767.30 MW of active power losses. Fig. 2 shows the evolution of the maximum active/ reactive power mismatches and Fig. 3 shows the evolution of the objective function. 5.3. Comparative tests Modi®ed Newton's method (MN), proposed here, had its performance compared with the following approaches: active set and penalty (ASP) [4], PD [5], and primal-dual logarithmic-barrier (PDLB) [6], in their coupled versions. All three approaches are Newton-based methods. The databases of the system had the same con®guration in all tests. The stop criteria took into consideration the convergence characteristics of each approach, and guaranteed a convergence within the region of feasibility and a maximum error for the equality constraints of less than 10 24. Further details are found in Refs. [4±6]. Table 4 shows the ®nal results of convergence for each approach on 810-bus system. The table presents the actives losses in the transmission, reactive power generation, the total number of iterations for convergence and the CPU time.
819
The results achieved in the tests show the viability of the use of the augmented Lagrangian in association with the modi®ed Newton method. 6. Conclusions ORD is a large scale nonconvex nonlinear programming problem with nonlinear constraints. The paper presents a new approach to the solution of this problem. The augmented Lagrangian corresponding to the original problem is solved by the modi®ed Newton method. The dif®culty in identifying the binding constraint set is removed by introduction of dual variables and quadratic penalty terms into the augmented Lagrangian. The method can start from infeasible points and verify the increase of the penalty factors, thus avoiding the ill-conditioned Lagrangian matrix. The tests demonstrated the effectiveness of the proposed approach. The hardest task in this approach is to ®nd the initial values of the penalty factors for the voltage magnitude. Acknowledgements This project was supported by FAPESP, FundacËaÄo de Amparo aÁ Pesquisa do Estado de SaÄo Paulo. References [1] Carpentier J. Contribution to the economic dispatch problem. Bull Soc Fr Elect Ser B3 1962;8:431±47. [2] Dommel HW, Tinney WF. Optimal power ¯ow solutions. IEEE Trans Power Apparatus Syst 1968;87:1866±76. [3] Sasson AM, Vitoria F, Aboytes F. Optimal load ¯ow solution using the Hessian matrix. IEEE Trans Power Apparatus Syst 1973;92:31± 41. [4] Sun DI, Ashley B, Brewer B, Hughes BA, Tinney WF. Optimal power ¯ow by Newton approach. IEEE Trans Power Apparatus Syst 1984;103:2864±75. [5] Santos A, da Costa GRM. Optimal power ¯ow solution by Newton's method applied to an augmented Lagrangian function. IEE Proc Gener Transm Distrib 1995;142(1):33±6. [6] Granville S. Optimal reactive dispatch through interior point method. IEEE Trans Power Syst 1994;9(1):136±46. [7] Costa GRM, Costa CEU. Improved Newton method for optimal power ¯ow problem. Electrical Power Energy Syst 2000;22:459±62. [8] Huneault M, Galiana FD. A survey of the optimal power ¯ow literature. IEEE Trans Power Syst 1991;6:762±70. [9] Hestenes MR. Multiplier and gradient methods. JOTA 1969;4:303± 20. [10] Luenberger DG. Linear and nonlinear programming. 2nd ed. Reading, MA: Addison-Wesley, 1984. p.436.
www.elsevier.com/locate/ijepes
Modi®ed Newton method for reactive dispatching G.R.M. da Costa Department of Electrical Engineering, Escola de Engenharia de SaÄo Carlos, USP, 13560-970 SaÄo Carlos SP, Brazil Received 22 December 2000; received in revised form 11 October 2001; accepted 13 December 2001
Abstract This paper describes a new approach to the optimal reactive dispatch problem, based on an augmented Lagrangian function of the original problem. The Karush±Kuhn±Tucker (KKT) optimality conditions are solved by the modi®ed Newton method. The second-order information in the original system of equations is approximated and the ®rst-order information is kept intact The proposed method requires less computer memory than those algorithms currently available. The effectiveness of the proposed approach has been examined by solving the IEEE 30bus and the BRAZILIAN 810-bus systems. q 2002 Elsevier Science Ltd. All rights reserved. Keywords: Power systems; Optimal power ¯ow; Optimization; Multiplier method; Augmented Lagrangian
1. Introduction The optimal reactive dispatch (ORD) is a particular case of the optimal power ¯ow (OPF) in which active controls are ®xed and control optimization only relates to reactive power, e.g. voltage magnitude for generators, synchronous condensers and static VAr, taps and switchable banks. The OPF problem was proposed by Carpentier in the early 1960s based on the economic dispatch problem [1]. Since then, many papers have been written in the attempt to solve it. Gradient techniques were the ®rst approach applied to the problem [2]. Slow convergence limited their use. Sasson et al. [3] proposed a method in which the objective function is penalized with both equality and inequality constraints. The Hessian matrix was used for the ®rst time. Sun et al. [4] introduced into OPF research a combination of the Newton's method with Lagrange multipliers and penalty functions. The major dif®culty turned out to be the identi®cation of binding constraints. The Dual-Newton approach [5] handles equality constraints with Lagrange multipliers and inequality constraints with both penalty factors and Lagrange multipliers. The binding constraints need not be identi®ed. In Ref. [6] a primal±dual (PD) logarithmic barrier algorithm is applied directly to a nonlinear OPF problem and the Karush±Kuhn±Tucker (KKT) conditions are solved by Newton's method. A crucial step in the algorithm is the choice of the barrier parameter. Costa et al. [7] showed a new approach that improves the performance of the Newton's method proposed by Sun et al. [4]. Up to the E-mail address: [email protected] (G.R.M. da Costa).
present there is no single robust, reliable and fast approach that ful®lls the requirements of the EMS. Huneault and Galiana [8] have made an extensive study of the approaches used to solve the OPF problem. It is evident that all these solutions have followed close behind advances in numerical optimization, and where possible were computed by techniques already well established, which use the incidence matrix structure. This structure has the advantage that computer packages can be applied to it. In this paper, we propose a new approach to the solution of the ORD problem. An augmented Lagrangian is associated with the original problem and the KKT optimality conditions are solved by the modi®ed Newton method. The constraints include the reactive power limits of generators, the limits on the voltage at the load buses, and the limits of the control variables, i.e. the transformer tap positions, generator terminal voltages, switchable reactive power sources and limits on transmission line ¯ow. All variables are treated as continuous. The paper is organized as follows: ®rst, the problem of ORD is shown. Next, the modi®ed Newton and the augmented Lagrangian approaches are discussed. The implementation of this new approach is described. The results obtained testing on two systems (30 and 810 buses) are reported. Finally, some conclusions are drawn.
2. Description of the ORD problem Minimization of real power losses in a system forms the basis of the reactive power optimization problem. The ORD
0142-0615/02/$ - see front matter q 2002 Elsevier Science Ltd. All rights reserved. PII: S 0142-061 5(02)00013-3
816
G.R.M. da Costa / Electrical Power and Energy Systems 24 (2002) 815±819
problem can be presented as:
constraints of the MORD model is de®ned by
minimize f x
L f x 1
i1
such that gi x 0 i 1; 2; ¼; m , n hj x # 0 j 1; 2; ¼; p
m X
(1)
where x [ R n is the vector of state variable, f(x) the real power losses of the transmission, g x 0 the power ¯ow equations, and h x # 0 is the limits on state variable and functional constraints of the power system. The state variable vector x represents the voltage magnitude, angles and LTCs taps. The objective function, f(x), represents the active power losses in the transmission. This function is nonseparable and permits no simpli®cations. The equality constraints, g x 0; represent the power ¯ow equations for scheduled load and generation. The inequality constraints, h x # 0; represent the state variables such as voltage magnitude limits and functional constraints of the power ¯ow, e.g. limit of active and reactive power ¯ows in the transmission lines and transformers, limits of reactive power injections for reactive control buses and active power injection for the slack bus. This is a typical nonlinear and nonconvex problem.
1
li gi x
p X v mj hj x 1 ej 1 hj x 1 ej 2 2 j1
where l i, m j are the Lagrange multipliers (related to the equality and inequality constraints, respectively). The minimum of L with respect to the slack variables (ej) is calculated from the following stationary conditions: 2L 0 2ej
4
that is:
mj 1 v h x 1 ej 0 j 1; ¼; p These conditions are satis®ed by mj ej max 0; 2 2 hj x v
p v X h x 1 ej 2 2 j1
such that gi x 0
6
MORD : p v X r 2 x 2 j1 j
7
such that gi x 0
When the problem (1) is analyzed in terms of its ordinary Lagrangian function, there is usually a duality gap, unless the objective and constraint functions are convex. The gap can be removed by passing to an augmented Lagrangian function that involves quadratic penalty-like terms. The modi®ed dual problem then consists of maximizing a concave function of the Lagrange multipliers and an additional variable, which is a penalty parameter. The problem becomes a sequential unconstrained minimization problem solved by Newton-based methods. The following modi®ed optimal reactive dispatch (MORD) model, related with the ORD, is obtained then as: MORD : minimize f x 1
5
Finally, the expression above is substituted in the MORD model, giving the following:
minimize f x 1 3. The modi®ed Newton method applied to an augmented Lagrangian function
3
2
hj x 1 ej 0 ej $ 0 where ej is the slack variable related to inequality constraint, hj(x) and v is the penalty factor. The classical Lagrangian corresponding to the equality
rj x # 0 where
8 > > < hj x;
mj v rj x m m > j > : 2 ; if hj x # 2 j v v if hj x $ 2
8
The Lagrangian related to the MORD model, as formulated in Eq. (7), is given by: p m X X v li gi x 1 rj x mj 1 rj x 9 L f x 1 2 i1 j1 The Augmented Lagrangian relating to the MORD model is de®ned by Eq. (9). The solution of the MORD model will be sought by an operational strategy de®ned on conditions imposed on the function L. The process involves the determination of values for x and l , such that the conditions of optimality are satis®ed. Thus: 2L 0; 2x
2L 0 2l
10
This system of equations was solved using the modi®ed Newton method in which Dx and Dl were calculated.
G.R.M. da Costa / Electrical Power and Energy Systems 24 (2002) 815±819
Thus, the Eq. (10) can be represented as: "
B
JT
J
0
#"
Dx Dl
#
" 2
7xL
slowly, to in®nity:
# 11
7 lL
where B is an approximation to the Hessian matrix, i.e. B is a diagonal matrix of second partial derivatives of L with respect to x. Each element of B is a linear combination of second partial derivatives of power ¯ow equations, its constraints and objective function. J is a matrix of ®rst partial derivatives of power ¯ow equations with respect to x. The second-order information in the original system of equations is approximated and the ®rst-order information is kept intact. There are various ways to choose the approximation B [10]. Here we assume B equal to the diagonal of the Hessian matrix, updated to each iteration. One advantage of this method is that B can be taken to be positive de®nite even when the Hessian matrix, is not as given in Ref. [10]. For conciseness, Eq. (11) can be written as: WB Dz 2b
12
where WB is the Lagrangian matrix, Dz the vector of corrections formed by Dx and Dl , and b is the gradient of L with respect to x and l . 3.1. Updating the multipliers and penalty factors The gradient of the augmented Lagrangian at the minimum point x p is: 2L 2f xp 2g xp 1l 2x 2x 2x 8 2h xp > > < m 1 vh xp 2x 1 > > :0
m v m p if h x # 2 v
if h xp $ 2
mk11 j
> > > :0
if hj x $ 2
mkj v
mkj if hj x # 2 v
vk11 bvk
15
where b is the penalty increasing factor (b . 1). The ill-conditioning usually associated with the penalty function approach can be controlled. The method treats the equality and inequality constraints in a uni®ed manner, and avoids a search for binding constraints. 3.2. Algorithm The algorithm proposed to solve the ORD problem is an iterative process consisting of the following steps: (i) Make starting estimates for z x; l and m x can be the same as the initial values for a power ¯ow l m 0 or any reasonable guess (ii) Evaluate b as a function of z (iii) If KKT conditions are satis®ed the problem is solved, otherwise (iv) Update m (Eq.(14)) and v (Eq.(15)) (v) Evaluate the matrix WB as a function of z (vi) Solve the system WB Dz 2b for Dz (vii) Update z by Dz (viii) Return to step (ii). The initial values allocated to the variables x can be in the infeasible region of the problem or can be a power ¯ow solution. In the evaluation of matrix WB in step (v), the updated values of m and v must be used. For suf®ciently large v, WB is positive de®nite, thus a problem can be convexi®ed (at least locally).
4. Implementation 13
The terms m 1 vh x can be seen as the Lagrange multipliers that satisfy the optimality condition of the problem (1), [9]. Therefore, a good approximation to the unknown m is: 8 > k k11 > > < mj 1 vhj x ;
817
Most of the work in the algorithm is in the solution of system (11). The Lagrangian matrix, WB, that results from the linear approximation of the KKT conditions, has a structure that facilitates the application of sparsity of techniques. The three-bus system of Fig. 1 is used as an example to show the structure of WB (Eq. (16)). Each symbol indicates a nonzero element and its coordinates indicate a speci®c derivative. For example, B with coordinates u2 ; u2 is 22 L=2u22 and J with coordinates V2 ; lP3 is 22 L=2V2 2lP3 2DP3 =2V2 : Only the horizontal
14
The penalty parameter v may also be adjusted during the process. As in ordinary penalty function methods, the sequences of v increase towards a ®nite value, or tend,
Fig. 1. Three-bus electrical system.
818
G.R.M. da Costa / Electrical Power and Energy Systems 24 (2002) 815±819
Table 1 Speci®cations of tests cases
Table 3 Optimization summary for 30-bus system
Systems
Buses
Lines
Q-controls
LTC-controls
30-Bus 810-Bus
30 810
41 1340
5 131
4 205
coordinates are shown: v1 2 B 6 6 6 6 6 6 6 6 WB 6 6 6 6 6 6 6 6J 6 6 6 4
u2
J
u3
v2
v3
lP2
lP3
J B B B B J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
J
lQ2 3 J 7 J7 7 7 J7 7 7 7 J7 7 7 J7 7 7 7 7 7 7 7 5
Iteration
Active power loss (MW)
Mismatch uDPuMax (MW)
Mismatch uDQuMax (MVAr)
1 2 3 4 5 6 7 8 9 10 11
18.76 18.38 17.80 17.41 17.40 17.33 17.26 17.25 17.23 17.23 17.22
0.1763 0.0091 0.0220 0.0163 0.0083 0.0023 0.0001 0.0000 0.0000 0.0000 0.0000
0.1411 0.0956 0.0867 0.0225 0.0237 0.0381 0.0034 0.0024 0.0006 0.0002 0.0001
16
This matrix is sparse and symmetric. It needs to compute and store only half of LU factorization due to symmetry. The matrix structure is constant through iterations, the ordering and symbolic analysis are done only once to create a static data structure. Thus the numerical factorization is carried out ef®ciently at every iteration.
multipliers related to the equality and inequality constraints are, respectively, l 0 and m 0: Table 2 presents the sensitivity of the optimization process to the penalty factors, v, of the voltage magnitude. On the basis of these results, the initial penalties were de®ned as v 500 for voltage magnitude and ratio tap while the other inequality constraints were v 1:0: All penalty increasing factors were de®ned as b 1:2: The process converged in 11 iterations. CPU time for this run was 1.14 s. The amount of reactive power generation was 147.16 MVAr, with a total active power loss of 17.22 MW.
5. Test results Tests were done to verify the ef®ciency of the proposed approach. The algorithm was implemented in FORTRAN, using double precision arithmetic on a SUN SPARC 20 workstation, in the Power Systems Optimization Laboratory of EESC, USP. The cases studied were the minimization of active power losses in transmission in the IEEE 30-bus and the BRAZILIAN SOUTH/SOUTHEAST 810-bus systems. Table 1 lists the speci®cations of these test cases.
Fig. 2. Mismatches for the 810-bus system.
5.1. 30-Bus system This test was accomplished with the following initial conditions: Vk 1.0 pu (100 MVA base) and uk 0:0 for k 1; ¼; 30; and ti 1:0 for i 1; ¼; 4: The Lagrange Table 2 Tests for 30-bus system Tests
b
v
CPU time (s)
Iterations
1 2 3 4
1.2 1.1 1.2 1.1
100 300 500 1000
1.71 1.63 1.14 Diverges
19 16 11 Diverges
Fig. 3. Objective function for 810-bus system.
G.R.M. da Costa / Electrical Power and Energy Systems 24 (2002) 815±819 Table 4 Comparative tests on 810-bus system Algorithm
Losses (MW)
Reactive (MVAr)
Iterations
CPU Time (s)
ASP PDBL PD MN
1690.01 1702.17 1689.42 1693.03
2498.73 2510.11 2496.89 2499.70
19 39 28 46
217.27 270.58 259.95 204.31
The optimization process for this case is summarized in Table 3. 5.2. 810-Bus system This test was accomplished with the following initial conditions: x t; V; u unconverged network solution. The Lagrange multipliers for the equality and inequality constraints were, respectively, l 0 and m 0: The initial penalties were de®ned as v 2000 for voltage magnitude and ratio tap and the other inequality constraints were v 1:0: All penalty increasing factors were de®ned as b 1:04: The process converged in 46 iterations to a speci®ed tolerance j 10 24 : CPU time for this run was 204.31 s. The amount of reactive power generation was 2499.70 MVAr, with a total active power loss of 1693.03 MW. The load ¯ow does not keep the voltage magnitude within desirable limits on the load buses, with a total of 2839.69 MVAr of reactive power generation and 1767.30 MW of active power losses. Fig. 2 shows the evolution of the maximum active/ reactive power mismatches and Fig. 3 shows the evolution of the objective function. 5.3. Comparative tests Modi®ed Newton's method (MN), proposed here, had its performance compared with the following approaches: active set and penalty (ASP) [4], PD [5], and primal-dual logarithmic-barrier (PDLB) [6], in their coupled versions. All three approaches are Newton-based methods. The databases of the system had the same con®guration in all tests. The stop criteria took into consideration the convergence characteristics of each approach, and guaranteed a convergence within the region of feasibility and a maximum error for the equality constraints of less than 10 24. Further details are found in Refs. [4±6]. Table 4 shows the ®nal results of convergence for each approach on 810-bus system. The table presents the actives losses in the transmission, reactive power generation, the total number of iterations for convergence and the CPU time.
819
The results achieved in the tests show the viability of the use of the augmented Lagrangian in association with the modi®ed Newton method. 6. Conclusions ORD is a large scale nonconvex nonlinear programming problem with nonlinear constraints. The paper presents a new approach to the solution of this problem. The augmented Lagrangian corresponding to the original problem is solved by the modi®ed Newton method. The dif®culty in identifying the binding constraint set is removed by introduction of dual variables and quadratic penalty terms into the augmented Lagrangian. The method can start from infeasible points and verify the increase of the penalty factors, thus avoiding the ill-conditioned Lagrangian matrix. The tests demonstrated the effectiveness of the proposed approach. The hardest task in this approach is to ®nd the initial values of the penalty factors for the voltage magnitude. Acknowledgements This project was supported by FAPESP, FundacËaÄo de Amparo aÁ Pesquisa do Estado de SaÄo Paulo. References [1] Carpentier J. Contribution to the economic dispatch problem. Bull Soc Fr Elect Ser B3 1962;8:431±47. [2] Dommel HW, Tinney WF. Optimal power ¯ow solutions. IEEE Trans Power Apparatus Syst 1968;87:1866±76. [3] Sasson AM, Vitoria F, Aboytes F. Optimal load ¯ow solution using the Hessian matrix. IEEE Trans Power Apparatus Syst 1973;92:31± 41. [4] Sun DI, Ashley B, Brewer B, Hughes BA, Tinney WF. Optimal power ¯ow by Newton approach. IEEE Trans Power Apparatus Syst 1984;103:2864±75. [5] Santos A, da Costa GRM. Optimal power ¯ow solution by Newton's method applied to an augmented Lagrangian function. IEE Proc Gener Transm Distrib 1995;142(1):33±6. [6] Granville S. Optimal reactive dispatch through interior point method. IEEE Trans Power Syst 1994;9(1):136±46. [7] Costa GRM, Costa CEU. Improved Newton method for optimal power ¯ow problem. Electrical Power Energy Syst 2000;22:459±62. [8] Huneault M, Galiana FD. A survey of the optimal power ¯ow literature. IEEE Trans Power Syst 1991;6:762±70. [9] Hestenes MR. Multiplier and gradient methods. JOTA 1969;4:303± 20. [10] Luenberger DG. Linear and nonlinear programming. 2nd ed. Reading, MA: Addison-Wesley, 1984. p.436.