OPTflAAL LOAD )FLOW SOLUTION IUSING TEE HESSIAN MATRIX

A.M. sson F. Viloria F. Aboytes

Department of Electrical Engineering
Instituto Tecnologico de Monterrey
Monterrey, N.L., MEXICO

grams. This, of course, does not imply that only work

along this line should be pursued, as pioneering work
The rapid convergence that Newton's method pos- done now may result in revolutionary ideas later when
sesses, by use of the Jacobian matrix, has led to an fully developed.
investigation of using a higher order matrix, the
Hessian, for an even faster convergence. It turns out Some time ago, one of the authors investigated
that this approach unifies the fields of nonlinear the of nonlinear programming minimization techni-
use
programming methods and Newton based methods. The ques to solve the load flow and optimal load flow prob
load flow problem can be defined as the solution of a lems.9"0 Compared with the use of Newton's aproach
system of simultaneous equations fi(x)=O, i=l, n. ..., on the same problems, it seemed that the two methods
It can be shown the Newton's method proceeds in a had very little in common, Newton's one beingby far
direction that minimizes F= 5 fi(x)2. The Hessian superior in computing time, storage and capacity of
load flow also minimizes F by assuming that it is a solving large problems. Nevertheless, work was con-
quadratic function, such that the linearizations tinued on the nonlinear programming approach. One
become HAx = - g, where the Hessian H is the matrix significant advance was made by introducing decomposi-
of the second partials of F and the vector g is the tion techniques""12 that made possible the solution
gradient of F. The optimal load f'low problem can be of large load flow problems with the approach. The
formulated by including some additional terms in F so optimal load flow problem, however, was found more dif-
that a single algorithm, based on the Hessian, essen- ficult to formulate in terms of a decomposed system.
tially solves both the normal and the optimal load flow
problems. An interesting aspect of the method is that Recent work along the nonlinear prograrmming ap-
an existing Newton'ss program can be updated to a proach has achieved the following:
Hessian program quite simply. The H matrix is some-
what less sparse than the corresponding Jacobian but 1) large load flow and optimal load flow prob-
enough so that sparse techniques should be used. lems can now be solved efficiently with comparable
Furthermore, the Hessian can be completely obtained computing time and storage requirements as Newton's
from the Jacobian, thus avoiding extra explicit func- method.
tion evaluations in the program. The paper presents
enough details of the method for an implementation of 2) surprisingly enough, the new ideas can be
a computer program. Numerical examples are given and programmed as an extension of existing Newton's load
compared with Newtons method results. flow programs.
INTRODUCTION 3) as far as the optimal load flow problem is
concerned, the new method has less convergence dif-
The introduction of Newton's method" for load ficulties than the one based on Newton's method8.
flow analysis, together with the use of sparsity tech-
niques, has -led to many new developments in recent The object of the present papers is to present
years. Among these, the following stand out: the method and to discuss it in terms of the first
new
two points just mentioned. The third point is present-
1) the use of minimization techniques2 to con- ed in more detail in another paper13 which compares
trol the convergence of Netton's method. several optimal load flow methods.
2) the use of different load flow equations3'4'6 NONLINEAR PROG_AIZWNG APPROACH
that result in faster convergence or less computing
time. Load Flow Problem
3) the semi-optimization formulations8'7 that Considering the load flow equations to be
make use of Newton's method to achieve area control and
maximum interchange capabilities. fi(x) =0, i1=, ..., n, (1)
4) the formulatiori of an optimal load flowvprob- where x=x1,..., xn, the minimum of the following func-
lem which optimizes production costs (economic dis- tion
patching), transmission losses, reactive compe nsation n
or other functions, with complete control on all sys- F(x) = Ifi(x)2 (2)
tem variables, such as voltage magnitudes and real and
reactive generation. i=l

equals zero and describes the solution of (1). Previ-

All in all, an impresive amount of work with pos- ous work9 11 used the Fletcher-Potvell method to mi-
itive results has beent done around the central basis nimize F(x). Starting from an initial point x0, a
of Newton's method. As a result of this, many util- correction ax is obtained by
ities have now a Newton's load flow program for pro-
duction purposes, although few take advantage of its AIx = - G g(x) (3)
full capabilities by using it as an optimal load flow
method. Future developments in load flow or optimal where G is initially a unit matrix and g(x) is the
load flow methods that build on the foundation of gradient of F, defined as
Newton's method stand a greater chance of being uti-
lized, specially if the new ideas can be applied by g (x) = aF (x) (4)
ax
extending or altering existing Newton's computer pro-
Paper T72 198-5, recommended and approved by the Power System Engineering Committee of the IEEE Power Engineering Society for presentation at the IEEE
Winter Meeting, New York, N.Y., January 30-February 4, 1972. Manuscript submitted January 8, 1971; made available for printing November 23, 1971.
31

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.
 The Fletcher-Powell method continues by varying the although symrmetric, that must be stored and altered at
matrix G based on information on the new point just every Fletcher-Powell step.
obtained and the previous one. With the new G, (3) is
again used to obtain a new correction for the x varia- 2) above 100 variables, the Fletcher - Powell
bles. This process is continued until F(x) approaches method breaks down in the sense that the alterations
zero with some tolerance. The power of the method lies of G fail to produce Ax values that reduce F(x).
in the fact that G approaches the inverse Hessian as
the minimum is located. The Hessian H is defined as 3) for the load flow problem this last charac-
teristic can be by-passed through decomposition, but
H =2 f (x) i,j=l, ..., n (5)
not so for the optimal load flow problem.
2xiXj 4) the number of Fletcher-Powell steps and the
and so, time per step increases with problem size.
G = H 1 (6)
THE HESSIAN LOAD FLOW
at the minimum point.
It is noted that the above comments point to two
If F(x) were a quadratic function in x, the ideas: (1), the weakness of the Fletcher-Powell meth-
Hessian would be a constant matrix. The correction od for large problems and (2), its necessity to work
with a full inverse Hessian matrix. To eliminate these
A x = - H-1 g(x) (7) two adverse conditions, it was decided to study the
implication of evaluating the Hessian directly. It
would take an initial point to the minimum inone step. turns out that for power systems applications, the
For an n variable quadratic F(x), the Fletcher- Powell Hessian of F(x) is a sparse matrix. Equation (7) can
method reaches the minimum in n steps. The salient then be written as
characteristic of the Fletcher-Powell method lies in
the fact that the inverse Hessian need not be computed, H Ax = - g(x) (9)
as an initial unit G is sequentially altered as the
minimization proceeds and only at the minimum does G and solved for A x by Gaussian elimination using spar-
become H_-1 sity techniques.

Optimal Load Flow Problem A solution algorithm for the load flow problem
could be as follows:
The optimal load flow problem can be stated as 1) an initial point is chosen arbitrarilly.
follows:
Minimize C (x) 2) the gradient and the Hessian of F(x) are com-
puted.
subject to constraints
I 3) the corrections ax are calculated from (9)
_;}(x) 3 0 , J =, .,1 using a Gaussian elimination scheme taking sparsity
where the number oz equality constraints nust be less into account.
Lnan the number of x variables. There is no Limit on
thne nunuber of inequality constraints. To solve the 4) steps 2 and 3 are repeated until F(x) ap-
ootirnal load flow problem, the following function is proaches zero with a prescribed tolerance.
constructed, m
From a practical points of view, the algorithmn is
F(x) = C(x) + E r. f.(X)2 (8) changed somewhiat. First, as soon as a Ax vector is
available, an a is found such that the correction /\x
j=l is the optimum one in the sense that it produces the
The solution proceeds as follows: smallest F(x 1 a Ax). By doing this, the method con-
trols convergence so that a non-divergent process is
1 given initCial r. values, F(x) is minimized obtained. This is the same idea that has been ap-
using the Fletcher-Powell inethod. The minimuri, of F(x) plied to Newton's riethod2.
is not now zero as in the 1load flow prob'lem. The
inequality constraints that are not violated are ex- The second alteration of the above algorithm has
to do with the way the Hessian and the gradient are
cluded from (8). computed. Appendix I shows that both the Hessian and
2) the r. values of violated inequality con- the gradient of F(x) can be computed directly from
straints are increased so that they have a grea'er Jacobian terms. Thus,if the Jacobian matrix is avail-
weight Jin (8). The same is done for equality con- able, no new explicit function evaluations are neces-
straints that do not approach zero f-ast enougn. sary.

3) steps 1 and 2 The third alteration has to do with the sparsity

are sequentially repeated until properties of tne Hessian. For load flow problems,
all constraints are satisfied to some tolerance and these properties can be fully determined from the to-
F(x) can be minimized no further.
pology of the system. For the Jacobian, it is well
Discussion on the Nonlinear Programming Approach known that if node i and j are connected, the bf i J
and the 6fJJxi are non-zero. For the Hessian, the
Experience obtained by using the nonlinear pro-
b2F/)xi,x. and the )2F/-3xj?xi are also non-zero.
However, it node i and j are not connected directly,
gramming approach has shown that the methods, as de- but through a third common node k, the two correspond-
scribed above, have some practical and theoretical ing Jacobian terms are zero but not so the two Hessian
deficiencies:* terms. Thus, the Hessian turns out to be less sparse
than the Jacobian, but sparse enough for sparsity tech-
1) the inverse Hessian matrix is a full matrix, niques to be exploited, for power systems applications.

32

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.
 Putting all the above comments together results sents the gradient and Hessian expressions.
in what is called in this paper the Hessian Load Flow
Method: In terms of power system variables, the problem
can be stated as follows:
1) an admittance matrix is computed and stored
in compact form. Minimize C(V,8), which is the cost of fuel or trans-
mission losses or any other function of system volt-
2) a connection table is formed which describes ages V and angles 6,
the non-zero element positions in the Hessian by ana-
lyzing the admittance matrix. subject to constraints,

3) an ordering scheme is used to optimally order P(V,I) Pe = f ixed real and reactive powers
the Hessian rows according to the information provid- at load busses.
ed in step 2. Q(V,6) =Qe
4) the Jacobian is computed and stored in com- Pgl P(V,6) <Pgu real and reactive generation
pact form. Actually, only certain off - diagonal between lower and upper limits.
Jacobian terms need be computed as the rest are either Qgl sQ(V,8) -:Qgu
equal to, or are the negative of, those computed.
V1 6 V % VU voltage magnitudes at all
5)
the Hessian and the gradient are computed one busses within prescribed limits.
row at a time and the triangularization process is 5,2E !~v
performed in compact form as soon as a row is avail- 8 =0O one bus angle selected as ref-
able. Back substitutions produce corrections A x. erence.

6) an optimal correction is determined by The load powers are thus equality constraints
finding an optimal a. while the rest are usually inequality constraints. It
may be that a particular case to be solved may have
7) steps 4 to 6 are repeated until F(x) is some fixed generation or voltage magnitudes. The meth-
smaller than some prescribed tolerance or until each od, as had already been explained, consists of taking
fi(x) is close enough to zero. all the constraints as the f. (x) functions of (8). A
lower and upper limit inequality are really two ine-
The reader will have detected by now that there quality constraints. The objective is to minimize
are points in common with a Newton's program: the ad- F(x), Equation (3), until all m functions fi(x) are
mittance matrix calculation, the ordering scheme and satisfied. Typical f. (x) functions are:
the Gaussian factorization on rows. In fact, a
Newton's program can be altered to work on the Hessian fp (V,8) = P (V 8) - PJ = 0 equality constraints for
algorithm. Step 2 would be an extra subroutine wihich a load at node j
prepares system configuration information so that the fqj (V,6) = Q (V, 8) -
Qj = O
same ordering subroutine now being used to order the
Jacobian would order the Hessian. Step 4 requires an
extra subroutine to store the Jacobian. The actual fpi (V,I) = P. (V,)-PglJI o lower inequality con-
straints for generation
calculation logic for the non-zero Jacobian terms is f qj (V I ) = Q (V,6)-Q gi
I j -;! at node j
the same as with Newton 's program. Finally, the sec-
tion of Newton's program that taIculates the Jacobian
by rows, operates on the rows and performs the back
f vj(V,8) =V. - V 1_0 lower inequality con-
straint for bus j voltage
substitutions, can be used here, provided a Hessian
row is computed instead of a Jacobian row. In power system terms,
N
The reader is referred to Appendix I for all the
expressions in terms of voltages, admittances and real
and reactive powers. It will only be mentioned here
F(V,8) = C(V,6) l (r f2
Pj +rqj f2 qjv
pjq +rf2 j ) (12)
j=
-

that the fi(x) functions are the real and reactive

power equations and the x variables are the bus volt- where N is the number of system nodes. It is impor-
ages and angles normally used by Newton's method. tant to stress that whenever an inequality constraint
Thus, the load flow equations are of the type, is within limits it is excluded from Equation (12).

fp(VM6) = P(V,8) - P required = 0 (10) To find the optimal solution the following proce-
dure is recommended:
f (V,6) = ((V,8) - Q required = 0 (il)
0) solve a load flow using the Hessian to get
where (10) is chosen for each generator bus, (10) and initial values for V and 6 variables.
(11) for each load bus and none for the slack bus. The
variables are 8 for aenerator buses and V and 8 for 1) chose all r values equal to one, initially.
load busses. The function to be minimized is made up
of the chosen equations, 2) perLormr two iterations of the Hessian Equation
(9).
F(V,8) = >fpi(V,S)2 + JiLqi(V,8)
ieG,L ieL examine all
3) constraints and multiply the r
terms corresponding to violated constraints by a
THE HESSIAN OPTflAL LOAD FLOW factor of 5 or 10.

For the Hessian optimal load flow, the problem is perform at most two more iterations of the
4)
stated in the form given by equation (8). Instead of Hessian. It mig'ht be that the /V and L8 corrections
using the Fletcher-Powell method to minimize F(X), the are too small in the first iteration and there is no
Hessian is evaluated and solved in the sanme way as need for the ne xt one.
described for the load flow problem. Appendix I pre-
33

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.
 5) if all constraints are within tolerance the Table II. Max. fi per Iteration for the Load Flow
optimum solution has been found. Otherwise, return to Problem.
step 3.
From the above presentation, it can be seen that
Max. fi - 14 bus Max. fi -30 bus
Itera- problem problem
the Hessian optimal load flow is just an extension of
the Hessian load flow. The F function is a bit dif- tion
Newton Hessian Newton Hessian
ferent now, and thus the expressions of the gradient
and the Hessian are altered by some additional terms.
Appendix I presents all of the terms of the Hessian. 0 .922 .922 .9345 .9345
The pattern of the Hessian in the optimal problem is 1 .122 .4649 .0631 .6678
similar to the one in the normal load flow problem, 2 .0009 .0408 .0005 .0847
and thus, all that was described in relation to ex- 3 lO .0002 10 .00l
tending a Newton's program for the normal load flow is 4 - .00004
also valid for the optimal problem. The word similar
has been used because in the optimal problem, genera- On first thoughts it might seem that the above
tor voltages may be variables and thus there will be results are contradictory. That is, a second-order
extra terms in the gradient and Hessian, although the method, as the Hessian, should be expected to converge
actual pattern is completely determined from the to- faster than a first-order one, as Newton's Jacobian.
pology in the same manner as in the load flow problem. However, it must be realized, that Newton's method in-
volves a linearization of the individual fi functions,
As far as storage is concerned, the optimal while the Hessian involves a linearization of the de-
Hessian method requires about the same amount as the rivatives of F =
optimal method based on Newton's method8. The Hessian
requires the storage of a triangularized nxn matrix
plus some terms of the Jacobian. The Newton's approach More rigorously, consider a general function ex-
requires the storage of a triangularized nxn Jacobian panded around an initial point xo
matrix plus a second triangularized nxn matrix to ob-
tain repeat solutions for the Lagrangian multipliers.
F(x) = F(x ) + g(xo) A x 4!bxt H(xo)ZA x +... (13)
As the Hessian is somewhat less sparse than the
Jacobian, the two methods require about the same order g(x) =)F =
g(xo) + H(xo)Lx +... (14)
c)x
of storage. In the Hessian approach, there is no need
for the Lagrangian multipliers calculations plus the If higher order terms are neglected, and if g (x)
alteration of the so-called control variables, which is set to zero, to describe the minimum of F(x),
is probably the weakest point of the Newton based op-
timal method. For a more detailed comparison of these 0 =
g(xo) + H(xo)ax (15)
two methods, reference is made to another paper.13
which can be written as
ANALYSIS BASED ON NUMERICAL EXAMPLES
H(xo)\ x = - g(x ) (16)
A bus system taken from chapter 10, refer-
5
ence 14, and the IEEE (AEP) 14 bus and 30 bus test This is Equation (9), used by the Hessian. In (13),
systems were chosen as examples as their data are H is the Hessian, the matrix of the second partials.
widely available. The above analysis shows that Equation (14) is the one
being linearized. This substantiates the statement
Table I. Convergence of the Hessian and Newton's Load that the Hessian linearizes the derivatives of F. If
Flow Methods. the individual fi are of power k, then F is of power
2k and the derivatives, of power 2k-1. If in the load
Iterations Optimum a/iteration flow problem, k=2, Newton's method linearizes
System a square function while the Hessian works on a cube
Newton Hessian 1 2 3 4 function. This accounts for the results of Table II.
Another item to bear in mind in a comparison, is
5 3 3 1.0117 1.1289 L.OC39 that Newton's method, as it is normally programmed,
14 3 3 .6836 .9%1 1.0273 requires an extra iteration to find out that it has
30 3 4 .4414 .918 l.1172 1.0078 reached the required precission. The Hessian does not
require this as the fi functions must be evaluated at
the beginning of each iteration.
Table I presents numerical results for the load The Hessian Optimal Load Flow
flow problem. It can be seen that both methods con-
verge in the same number of iterations and that the As examples, the Hessian optimal load flow method
corrections are substantial in the Hessian. Actually, was used to solve the same three systems used as exam-
Newton's method converges a bit more rapidly than the ples of the Hessian load flow. Tables III, IV and V
Hessian. This was evidenced by examining the greatest present the final results plus the limits used and the
fi for both methods. All fi should be zero at the so- cost function for each case.
lution according to the load flow problem formulation
given by Equation (1). The maximum fi value at each The convergence characteristics are summarized
iteration is a measure of the quality of convergence in Table VI. It can be seen that a very rapid con-
and, of course, of precission of results. Table II vergence is obtained in all cases with the maximum
presents this comparison. constraint violation under control at all times.

34

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.
Table III. Hessian Optimization Results, 5-Bus Case. Table VI. Hessian Optimization Convergence Character-
istics.
NVmin V Vmax 8 min IPmax min Q Qmax
1 1.02 1.02 1.02 0 0.3 .954 1.2 0.0 .27 .6 Iteration 1 Iteration 2 Iteration 3
2 0.9 .955 1.05 -6.32 .6 .3
3 1.04 1.04 1.04 -1.98 0.3 .697 1.2 0.0 .527 .6 tSem Maximum Maximum -Maximum
4 0.9 .923 1.05 -9.40 .4 .1 F Viol. F Viol. F Viol.
5 0.9 .993 1.05 -4.08 .6 .2
5 771.12 1.OX10-3 763.05 1.35X10-3 760.09 1.5X10-3
C = 100
200
p + 3.51P +44.4+ 21
1 p+ 3.89P3+
* 3
40.6 14 1154.9 4.88X16 1146.8 2.OX103 1142.7 2.29X10
30 1260.9 5.6X104 1245.8 1.99X10 1245.6 2.OX10-3
Table IV. Hessian Optimization Results, 14-Bus Case.

_ Vmin V vmax 6 Pmin _ max Qmin Q Qnmax

Sys- Iteration 4 Iteration 5 Iteration 6
1 1.06 1.06 1.06 0. 0.5 1.58 2.0 -.2 -.c03 1.0 tem Maximum Maximum Maximum
2 1.0451Q045 1.045 -3.2 0.2 .665 1.0 -.4 .204 0.5 F Viol. F Viol. F Viol.
3 1.01 1.01 1.01 -10.4 .942 -.4 -.02 O.6
4 0.9 L029 1.1 -7.7 .478 -.04 5 760.9 1.3X10-3 760.8 1.2X10-3
5 0.9 1.037 1.1 -6.3 .076 .016 14 1137.0 l.lX103 1136.4 3.5X10 3 1136.4 1.X10-4
6 1.07 1.07 1.07 8.4 .2 .433 1.0 -.06 .412 .45 30 1242.4 1.6X10- 1242.5 1.5X10-3
7 0.9 1.04 1.1 -9.8 0.0 0.0
8 1.09 1.09 1.09 -9.8 0.0 -.06 .l11 .45
9 0.9 L 023 1.1 10.9 .295 .166
10 0.9 1.02 1.1 -10. 8 .09 .058
11 0.9 .04 1.1 -9.7 .035 .018
12 0.9 .05 1.1 -9.4 .061 .016 The initial value used for the r terms should be
13 0.9 .0451.1 .6 .135 .058 chosen so that there is a balance between the initial
14 0.9 L014 1.1 11.31 .149 .05 values of the partials of C and those of the con-
straint term in Eq. (10). This is very important as
C 200+2 .45P1 2005L+3.51P2+44.4+21 +3.89P6+40*6 poor convergence results if these two quantities are
different by several orders of magnitudes. Care must
Table V. Hessian Optimization Results, 30-Bus Case also be exercised in not increasing the r terms by
too big a factor during the process. By doing this
convergence might terminate prematurely at a sub-
N Vmin V max 6 Pmin P Pmax mn optimal point. Factors of 5 to 10 were found accept-
able.
1 1.06 1.06 1.06 0. 0.5 1.7 3.0 -.6 .013 .6
2 1.045 1.045 1.045 3.3 0.2 .738 1.5 -.6 .397 .6 CONCLUSIONS
3 0.9 1.0271.1 -5.5 .024 .012
4 O.9 10191.1 6.6 .076 .016 The paper has presented the underlying theory be-
5 1.01 1.01 1.01 -11.5 .942 -.6 .381 .6 hind the Hessian method of solving the load flow and
6 0.9 1.013 1.1 -7.8 .0 0.0 optimal load flow problems. Also, results of applying
7 0.9 1.003 .1 9.8 .228 .109 the method to actual systems have been included. All
8 1.01 1.01 1.01 8.4 .3 -.6 .313 .6 this shows that the method has many attractive charac-
9 0.9 1.025 L1 7.5 .0 0.0 teristics:
10 0.9 1.0051.1 >10.4 .058 10.2
11 1.082 1.0821.082 2.1 0.2 .497 1.5 -.6 i.3172 .6 1) comparable convergence and storage properties
12 0.9 1026 1.1 -11.4 .112 1.075 to Newton's method for the load flow problem. However,
13 1.071 1.071 L 071 11.4 0.0 -.6 .344 .6 the method cannot be justified for this application
14 0.9 L009 1.1 12.2 .062 .016 only.
15 0.9 1.003 1.1 12.1 .082 .025
16 0.9 1.01 L1 11.3 .035 .018 2) with the Hessian method, the optimal load
17 0.9 1.002 1.1 -Q9 .09 .058 flow problem becomes a problem of the same type as
18 0.9 .992 1.1 -12.2 .032 .009 the normal load flow problem. Convergence for the
19 0.9 .988 1.1 12.1 .095 .034 optimal problem is thus very fast.
20 0.9 .992 1.1 117 .022 .007
21 0.9 .992 1.1 -11.0 .175 .112 3) an existing Newton's load flow program can be
22 0.9 .992 1.1 11.0 0.0 0.0 extended easily to become a Hessian load floxi and op-
23 0.9 .988 1.1 12.1 .032 .016 timal load flow program.
24 0.9 .977 1.1 -11.8 .087 o.67
25 0.9 .975 1.1 -12.0 0.0 0.0 4) the method can be applied efficiently to
26 0.9 .957 1.1 -12.5 .035 .023 large systems.
27 0.9 .983 1.1 11.9 0.0 0.0
28 0.9 1.008 1.1 -8.3 0.0 0.0 These characteristics imply that while the method
29 0.9 .963 1.1 13.2 .024 .009 cannot be justified only as another load flow method,
30 0.9 .95 1.1 -14.2 .106 .019 it is extremely powerfull as an optimal load flow
method that can be also used to solve the normal load
flow as a by-product. In fact, the authors venture in
C 200P+2.45P1+l05+Pe+3.51P2+44.4+2fK+3.89E1+40.6 saying that it is the most attractive optimal load
flow method known to them. A comparison of the
Hessian method with other optimization methods13 tends
to back this strong statement.

35

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.
 APPENDIX I RPj(APjNji - N. N..) + RQ(!Q H. . - H. . L.) +
i i ii .1j31 3i 331 J133
This Appendix presents the expressions for the
gradient and Hessian of F(V,8). The notation is as RPi(L PiNij -NiiNij)+RQi (QHij -Hij Lii ))+V v
follows:
N N
Jacobian relation assummed to be known, (5) V = -2(> RPk(-PkHki-NkiHki) RQk
1j ~k=l k=l
H N; ' | i=j
LJ L IV/V Q
(AQkJki LkiJki)) + Vi 2)vi
a
RP = r factor related to an f constraint.
RQ = r factor related to an constpraint. N N
RV = r factor related to an f constraint. a2F
The function to be minimized is C(V,8), a general
(6) Vi i2T
BV = - 2(- l
k=l RPkNkiHkj+ k=l RQkHkiNk.
> +

k?i k#i
function which could be the total cost of production, i#j k#j k#j
the interchange between two areas, or any other power
system function to be optimized. In the equations that
follow, no expressions are given for the derivative of
C(V,8) as they depend on the particular problem being RPi(Lb PiHij-NiiHij) + RQi(LiiNij -AQiNij) - RP
solved.
(AP H..+N..H..)+RQ.(&Q N. R-H J .))+V.
There are two expressions for the gradient,
N N

av i ( 7
k=l
RPk LPkNki+
k=l
QkLki+ TI RQkl REFERENCES
RViV iVi) + Vi C (1) W.F. Tinney and C.E. Hart, "Power flow solution
av. by Newton's method," IEEE Trans. PAS-86, pp. 1449-
N N 1460, 1967.
RQk AQkJki) +;-.
)F
(2) =- 2 ( Z RPk APkHki+ (2) A.M. Sasson, C. Trevifio and F. Aboytes, "Improved
1)8 k=l k=l Newton's load flow through a minimization tech-
nique," IEEE Trans. PAS-90, pp. 1974 - 1981,
Hessian 1971.
There are six expressions for the Hessian, (3) H.W. Dommel, W.F. Tinney and W.L. Powell, "Fur-
ther developments in Newton's method for power
~2aF\ N N system applications ," Paper No. 70-CP-161-PWR
presented at the IEEE Winter Power Meeting, 1970.
(1) = 2( LI RPk(PkNki+Hi)+ RQ
i=j k#i k#i (4) C. Trevifto, "Cases of difficult convergence in
load flow problems," Paper No. 71-CP- 62-PWR
)2c presented at the IEEE Winter Power Meeting, 1971.
+ RPi(Hii +APiJii) + RQi(J4i-&QiHii)) +
(5) Y.P. Dusonchet, S.N. Talukdar, H.E. Sinnott and
N N A.H. El-Abiad, "Load flows using a combination of
point Jacobi and Newton's method ," IEEE Trans.
(2) = -
(
3k=l
(LRP kHki
ki3 -Hk k=l RQkNkk
Nk PAS-90, pp. 941 - 949, 1971.
i j k-i k#i
k#j k#j (6) J.P. Britton, "Improved area interchange control
for Newton's method load flows," IEEE Trans. PAS-
+ RP.(\PiNji HjiHj.)RQi(L&QiH.i + j..Nj i)
- +
88, pp. 1577-1581, 1969.

+ RPi(&PiN.j H iiHj) + RQi(QHQiHi;+ Jii.N.

- (7) J.P. Britton, "Improved
through a more
load
general equation form,"
flow performance
IEEE
Trans. PAS - 90, pp. 109 - 116, 1971.
13 N N
(8) H.W. Dommel and W.F. Tinney, "Optimal power flow
= - 2(- RPk(N2 -&LP N + Q solutions," IEEE Trans. PAS-87, pp. 1866-1876,
i1jA.VO.
(3V
ki k ki k
1 3 ~k=l k=l 1968.
i=j k#i k#i
(9) A.M. Sasson, "Nonlinear programming solutions for
(AQkHi -H2i) -RP.(N2i.- &Pi(Nii+2ViGii))-RQ.(Li.- the load flow,minimum loss and economic dispatch-
ing problems," IEEE Trans. PAS-88, pp. 399-409,
1969.
AQi(L..-
1 11 2V21 B..))+RV.(V
N 1
-
iVi)V
i i ) +V.V.i
iOi (10) A.M. Sasson, "Combined use of the Powell and
)2F+
() i j 3V = 2(-
k RPkNkjNki
k RQkHkjHki Fletcher-Powell nonlinear programming load flow
1- 3 V k=l method," IEEE Trans. PAS-88, pp. 1530 - 1537, 1969.
kl
, . k7&i ksi
3.#J k#j k#j

36

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.
(11) A.M. Sasson, "Decomposition techniques applied to Thus it can be seen that for the system equality constraints the
the nonlinear programming load flow method," IEEE Lagrange multiplier will not necessarily be zero. A comparison with the
Trans. PAS-89, pp. 78-82, 1970. derivative of the authors' equation shows that:

(12) F. Aboytes and A.M. Sasson, "A power systems de- 2rj f (x) = 'j for j = 1, . . ., m
composition algorithm," 1971 PICA Conference
Proceedings, pp. 448 - 452, 1971.
For those inequality constraints which are not limiting Xj = 0, which
implies that rj must also be zero and of course this corresponds to leav-
(13) A.M. Sasson, F. Aboytes, R. Cardenas, F. Gomez, ing those constraints out as the authors indicate. Similarly, for those
and F. Viloria, "A comparison of power systems constraints which are limiting, X-j * 0, which implies that 2r fj(x) * 0,
static optimization techniques," 1971 PICA and since fj(x) is to approach zero, rj must become correspondingly very
Conference Proceedings, pp. 328 - 337, 1971. large if the two methods are to approach the same optimum. This may
cause numerical problems as indicated by the authors.
(14) W.D. Stevenson, Jr., Elements of Power System My point in drawing the above comparison is to show that the
Analvsis. 2nd. edition. New York, Mc-Graw Hill authors' method is similar in certain respects to classical Lagrange mul-
1962. tiplier techniques. The chief drawback to these methods is the way in
which inequality constraints are handled. The method chosen by the
authors consists of solving a succession of optimization problems, hav-
ing only equality constraints, by switching the inequality constraints in
and out, as they become violated or satisfied. In the classical Lagrange
multiplier approach the Kuhn-Tucker conditions indicate which con-
straints are to be imposed or ignored and a similar approach could be
used in the authors' program. However, the Kuhn-Tucker conditions,
when used as constraint switching rules, do not guarantee convergence
Discussion to an optimum. At each iteration in such a process there are a large
number of combinations of constraints which could be imposed or
ignored and no indication which ones may send other constraints into
B. F. Wollenberg (Leeds and Northrup Company, North Wales, Pa. violation if relaxed, the net results, especially for large systems, can be
19454): A good deal of insight into the authors' proposed optimal load an aimless search taking quite long or in fact never converging. It is for
flow solution can be obtained by a direct comparison with the classical this reason that optimization techniques such as linear programming
Lagrange multiplier technique. First, the author's equation (8) can be and quadratic programming have come into wide use. In LP and QP the
compared to the same problem formulated with Lagrange multipliers: rules for constraint switching are much better defined, the objective
function is always reduced, and the iterations always move from one
feasible point to another. Convergence to the optimum is virtually
Author's Optimal m guaranteed.
Load Flow Equation: F(x) = C(x) + Z rij(x)2 In looking at the authors' test cases it can be seen that none of the
j=l inequality constraints are limiting at the optimum - the resulting rapid
convergence is not surprising. Have they tried any cases with many con-
constraints limiting at the optimum and if not do they plan to? Without
Lagrange Multiplier m such tests the usefulness of the authors' method cannot be assessed.
Formulation: F(x) = C(x) + 1 Xj fj(x)
j=l REFERENCE
[1] Lasdon, L. S., Optimization of Large Scale Systems, MacMillan,
The optimum is achieved by setting the first derivatives of each equa- 1970, pg. 77.
tion to zero.

First Derivative of aF(x) =C(x) afj(x) A. Rothman (15, Rue Pierre Semard, Paris, France): The authors are to
Authors' Equation: ax axi + 2rifi(x) a x 0
be commended for having advanced with a very important step the load
' j=l
flow optimization solution. Dr. Sasson contribution in solving optimal
load flow problems via a nonlinear programming method should be par-
First Derivatives of aF(x) aC(x) fj() 0 ticularly emphasized. [ref. 9, 10, 11 and 13 of the paper] .
Lagrange Equation: axi ax),, . . ., n In eq. (3) of paper [9] Sasson has pointed out that the correction
x of the vector x (the so called activity vector) in an unconstrained non-
linear optimization problem, is computed by
aF(x) = =0 j= . m
Ax]=[-H 1] ax
The optimality conditions on the Lagrange equation are the well known
Kuhn-Tucker conditions:1
To minimize C(x), subject to f(x) > 0 j = 1,...,m where H = [af(x) / (axaxj) ] is the Hessian
Solve for: and g = [af(x) /ax] is the gradient column matrix.
m
VC(x°) + Z Xj (Vt(x') = 0 In their present paper emphasis is placed on a new sequence of
j=l computation in which both H and g are used for computing Ax consider-
ing at the same time instead of the unconstrained function f(x) a Kuhn-
and 0 if ((x)> O Tucker function F(x), adjusted in accordance with a previous method
used by Sasson in his application of the Fletcher Powell technique.
There are two additional important advantages of the method
or proposed:
1. It can be derived as an extension of Newton's programme i.e.
<°
S 0 if fj(x) = 0 by using the already computed Jacobian matrix.
2. It solves also the load flow problem as a by-product.
The problem as envisaged by the authors is general and mathemat-
where x°, X° are the values of x and X at the optimum. ically correct and complete that is it represents up to this day the most

Manuscript received February 17, 1972. Manuscript received February 22, 1972.

37

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.
adequate model of the optimal load flow. However I should like to H. Duran (University of the Andes, Bogota, Colombia): The authors
point out that in all practical applications consideration should be given have presented a new approach to solving the optimal power flow prob-
to the following conditions: lem. To obtain an initial point for the optimization process, a load flow
- the input data should not be too numerous and their accuracy is solved. In the paper, it is shown that the Hessian optimization method
should be compatible with the built in conditions of the problem. can also be used to solve the initial load flow. The authors further point
- the constraints introduced should be carefully balanced with the out that numerical experience has indicated that a Newton program
effective conditions of the power system so that they might lead to an would be faster than the Hessian for the load flow. This discussion pre-
exact interpretation; in other words constraints for instance related to sents a clear demonstration of these convergence behaviours for the
active and reactive branch loads, and which are presently in general simple case of solving one equation with one variable:
neglected, should play a major role as they represent important physical
limitations which when accounted might greatly improve the conver- f(x) = 0 (1)
gence of the solution.
I should like to ask the authors:
- the total number of programme instructions. Correction by Newton's method
- the total computer time for solving the 30-bus case.
- an appreciation of the time for solving a large network (say 400 xk+l = xk +Axn (2)
nodes).
where
Axn = - f (3)

S. T. Despotovic (Electrotechnical Research Institute, Belgrade, Yougo- Correction by the Hessian method
slavia): It is well known that the minimum of the quadratic form in x xh = xk+
xk+l h Axh
(4)
n where
F(x) = E f(x) Axh = - g'
i=l (5)
g= f.f' (6)
is the solution of the load flow equations
g, = (f)2 + ff' (7)
fj(x) =0, i 1. n
Substituting Eqns. (6) and (7) into (5),
x= XI, ...* Xn
Axh
Ah= - 1' ff (8)
Starting from initial point xo the authors solve the expression
Finally, substituting Eq. (3) into (8),
H-(x-xo) = -g(x)
Axh = Axh (l + ff)) (9)
for x-x by Gaussian elimination using spasity techniques. The solu- (f 1)2
tion x isihe unique minimum of the F(x) and at the same time the solu-
tion of the fi(x), i.e., of the load flow problem. The authors have given There are two cases to consider:
the name to the procedure of solving load flow problem the Hessian
load flow. Case]: f>0, f">0 or f<0, f"<0 (see Fig. 1)
For solving the optimal load problem the authors use the Hessian In this case, Eq. (9) shows that
load flow and Zangwill's method.
The Hessian method of solving the load flow problem does not lAxhl < lAxnl
appear attractive at all. But the combined Hessian and Zangwill's meth-
ods for solving optimal load flow problems is very attractive.
It would be, perhaps, better to solve a load flow using, instead of which clearly indicates that the Hessian converges slower.
the Hessian, Newton's method, in order to get initial values for the V
and 0 variables, and then to find the optimal solution by the Hessian Case 2: f > 0, r < 0 or f < 0, f" > 0 (see Fig. 2)
optimal load flow. The reasons are the following: In this case, Eq. (9) shows that
1. Newton's method is faster.
2. For the Hessian load flow, the Hessian matrix and the gradi- IAxhI > IAxnI
ent are of

Z
n

i=l
f(x)2

and for the Hessian optimal load flow of

f1
which also clearly indicates that the Hessian converges slower.

x x

k+l Xk x

n Fig. 1. Case 1
C(x) + fi(x)
i=l
f
In any case the Hessian load flow and Hessian optimal load flow
are superior to the methods described in References 9 and 10. The au-
thors should be commended for their efforts in developing methods for xk+l Xk x
solving power system problems. Fig. 2. Case 2
Manuscript February 22, 1972. Manuscript received February 22, 1972.
38

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.
0. A. Klitin (American Electric Power Service Corporation, New York, S. S. Sachdeva and R. Billinton (University of Saskatchewan, Saskatoon,
N.Y. 10004): I would like to congratulate the authors for a very well Sask., Canada): The authors are to be congratulated on an interesting
written and interesting paper for solving the optimal load flow problem. paper. We agree completely with their observation that it is much easier
These comments are specifically directed to the conventional load to modify an existing load flow method for optimal load flow solution
flow solution which must be obtained prior to starting the optimization than to evolve a completely new method and put it into practice. In
process. It is shown in Table I and in the paper that a Newton Load their earlier paper, the authors introduced a scaler multiplier oa to con-
Flow has better convergence characteristics than a Hessian load flow. trol the convergence of Newton's method. In the present paper, the
Nevertheless, it might be possible to structure the Hessian load flow to authors have introduced another step for evaluating the Hessian matrix
obtain Newton convergence characteristics as follows: and gradient from the Jacobian terms. The gradient and Hessian are
The Newton load flow solution is given symbolically by: then used in the following equation to determine the required change in
Ax.
Ax = J-1Az (1) HAx = -g (1)
J = Jacobian
Az = z - f(xk)
Equation (1) is obtained by expanding the function near the optimum
point using Taylor's theorem. The solution is obtained by Gaussian
AX = Xk+l - xk elimination techniques utilizing the sparsity of the Hessian matrix H.
The Hessian optimal load flow technique is a useful addition to the
literature because it utilizes part of the Newton load flow algorithm in a
Equation (1) can be rewritten as: non-linear programming area.
In the Fletcher Powell method, the correction vector Ax is given
by the following equation which is essentially the same as equation (1).
Ax = (JTJ)-l JTA (2)
Ax = -H-Wg (2)
Equation (2) is mathematically equivalent to (1). Numerically this
equivalence should hold assuming J is well conditioned. The resulting The Hessian and Fletcher Powell approaches differ only in the computa-
formulation results from expanding the authors' scalar function in a tional details required to solve equations (1) and (2). A comparison be-
slightly different fashion. tween the Hessian approach and the Fletcher Powell technique from
which the Hessian method is obtained, brings out some of the compu-
n tational requirements of the two approaches. The basic computational
F(x) = Z (zi-f,(x))2 (3) steps of the two methods are as follows.
i=l Comparing the main steps it may be seen that the authors have
used block CH in Figure 2 instead of block CF used in the Fletcher
Powell approach as shown in Figure 1. In block CF the matrix H is built
Where z. is the data for the (ith) load flow equation. Rewriting F(x) as a by the blocks CF2 and CF3 in which additional n X n matrices Ml and
dot proAuct of 2 vectors we get: Ml are evaluated n times to determine the optimal point for a quad-
ratic function. The number of iterations for the power system functions
are more than n. In the Hessian method CF is replaced by CH and H-'
F(x) = [z-f(x)]T [z-f(x)] (4) or Ax is evaluated at each step using a Gaussian elimination technique.
In the Fletcher Powell method, the function and gradients (Block BF)
are calculated using analytical gradients of the penalized function. In
Expanding f(x) in a Taylor series, truncating after the linear term we the Hessian approach the Jacobian matrix is used to evaluate these
have: terms (Block BH). Blocks AH and AF are virtually the same. It would
be of interest to know the computational effort required for each of
the two methods to solve the same problem. In particular we would also
F(x) = [Az-JAx]T[Az-JAx] (5) appreciate an indication of the requirements for the Hessian inverse
H-1 using Fletcher Powell and by Gaussian elimination techniques as-
F(x) = AzTAz - 2AxTJTAz + AxTJTJAx (6) suming that the two methods require almost the same effort to carry
out blocks A and B. A principal point is in regard to which of the two
methods requires the least computational time and money requirement.
The above quadratic form can be compared with the authors' expansion It would also be of interest, to know the maximum size of the problem
term by term. It can be shown that although jTj H, they are struc-
= solved using the Hessian approach.
turally identical. The minimum of the quadratic form (6) is: It appears that the main attraction in using a Gaussian elimination
technique to solve equation (2) for Ax is the sparsity of the matrix H. If
the sparsity [3] is judged as k/nn with k the number of non-zero ele-
ments in the n X n matrix, it would be of interest to know the sparsity
Ax = I(JTJ)-I (_ 2JTAz) (7) of H for the various systems studied. The authors comments on any
computational experiences regarding the round-off error associated with
the Gaussian elimination would also be appreciated.
The authors' g(x) is the gradient of eqn (4): The choice of the initial values of the penalty constants r depends
upon the violated constraints. As noted in the paper they do effect con-
siderably the convergence stability. Guide lines or suggestions for the
g(x) = -2JTAz (8) initial selection of the penalty constant r with regard to the violated
constraint would be of value. Non-linear programming is a powerful
analytical tool for solving many types of decision making problems.
Substituting (8) into (7) we have Variations of these techniques which attempt to solve a larger problem
can be very useful. In an attempt to solve larger optimal power system
operational problems efficiently the decomposition technique based on
the variable decoupling of the power system mathematical function
Ax -
= I(jTj)71 g(x)
2 have been proposed [1, 2]. It is believed that the application of these
techniques can make the suggested Hessian approach more efficient.

Thus, the initial load flow could be calculated assuming H = 2(JTJ) REFERENCES
without altering the logic to their present program.
I would appreciate the authors' comments with regard to this ap-
proach and possibly to expand Table I to include numerical results using [1] R. Billinton, S.S. Sachdev, "Real and Reactive Power Optimiza-
this approximation to the Hessian matrix. tion by Suboptimum Techniques", 71 CP 596-PWR, IEEE Summer
Power Meeting, July 197 1.

Manuscript received February 22, 1972. Manuscript received February 24, 1972.

39

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.
 AF
Initial vect( Dr To find the argument To find new gradient
and gradient vector x f(x ) } at x from analytical
x° aind go
= min f (x g ) expressions
-P ¢ ~~~n-cycles
_i.
t
I To find new i; T+lL.JTo find correction To find correct1ion
I H i+l=
I I4
H+M +M 2F matrix Ml: matrix M1:
1
I -H1Ag1Ag H1 i Ax'Ax' I
M = I
I 1 .tI
I Ag H1Ag1 Ax1 AgI1 I
I
F3 F F
- - - - - - -

CF
Fig. 1. Fletcher Powell Method

B AH
Initial vector x°-..3 To find the argument
vector xi l
+
f(x )= mnf(x
- aAx )

Using Gaussian elimina-

tion find correction
1vector Ax i :.H ii
Ax = -g i
CH
Fig. 2. Hessian Approach

[2] C.H. Jalissaint, N.V. Arvantidis and D.G. Luenberger, "Decompo- same as that of the Hessian approach. Another important point, is that
sition of Real and Reactive Power Flows: A Method Suited for On- the Fletcher-Powell method deals with the inverse Hessian matrix, usu-
line Application" IEEE Winter Power Meeting, 71 TP 113-PWR, ally a full matrix, while our approach exploits the sparsity of the Hes-
1971. sian matrix. This condition produces an excessive storage requirement
[3] E.C. Ogbuobiri, W.F. Tinney and J.W. Walker, "Sparsity-Directed with the Fletcher-Powell method as system size increases.
Decomposition for Gaussian Elimination on Matrices", IEEE Trans. The choice of the initial values of the penalty constants r is such
PAS-89, January 1970. that the two components of the penalized objective function, Eq. (8),
have the same order of magnitude. Initially, all r terms are given the
same value.
The sparsity of the Jacobian and Hessian matrics are given below
for various system sizes (% n K n X 100, where K is the number of
nonzero terms and n the system size):
Albert M. Sasson, Franklin Viloria, and Florencio Aboytes: We shall
first address ourselves to Mr. Klitin's discussion because his load flow
formulation improves greatly the performance of the first step of our SYSTEM JACOBIAN HESSIAN
method. With his idea, the iteration matrix for the load flow preserves 14 29% 57%
the sparsity of the Hessian while retaining the convergence charac- 30 13% 28%
teristics of Newton's method. We had justified a less efficient initial 57 6.6% 13.8%
Hessian load flow against using Newton's method, to avoid ordering 118 3.5% 9.2%
equations twice, once for the Jacobian sparsity and again for the Hes-
sian. Mr. Klitin's idea has been tested and reported elsewhere [ 15 ] .

Mr. Despotovic's suggestion of using Newton's method for the ini- As can be seen, although the Hessian matrix is approximately twice less
tial load flow are answered by our comments to Mr. Klitin's discussion. sparse than the Jacobian, the sparsity levels for larger systems are quite
Mr. Duran also addressed himself to the initial Hessian load flow low, warranting its exploitation by use of optimal ordering techniques.
and shows that, for the one-variable problem, Newton's method ap- Referring to Mr. Rothman's comments, the data for the optimal
proaches that solution faster than the Hessian. Again, with Mr. Klitin's load flow is much the same as that for an ordinary load flow, with the
idea, this unwanted situation is by-passed. exception of cost coefficients and limits on variables. Line flow con-
Referring to Messrs. Sachdeva and Billinton's discussion, we would straints can easily by taken into account in the formulation, but we have
like to recall again some of the differences between the Hessian ap- not included these constraints in our tests.
proach and the Fletcher-Powell method as stated in the paper. One main With respect to Mr. Wollenberg's discussion, we were aware of the
difference between the two approaches is that in Fletcher-Powell, the relationships he presents. In fact, in Reference 13 of the paper we dis-
inverse of the Hessian matrix is built iteratively, starting from a positive cussed the agreement of our approach with the Kuhn-Tucker conditions.
definite matrix, usually the unit matrix, while in our case, the Hessian These conditions are just that, and not a method of solving a problem.
matrix is directly obtained calculating the second partial derivatives. Thus the difficulties that Mr. Wollenberg states as to a guaranteed con-
Therefore, the correction vector Ax given by Fletcher-Powell is not the vergence to an optimum. However, we do not share his view that for
large systems, constraints will necessarily move aimlessly into, and out of,
Manuscript received April 26, 1972. a violation state. This depends on the method used to perform the mini-

40

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.
mizations and by the speed by which r values are increased. Assuming r are always present, we believe our results to be a valid indication of the
values are slowly increased, our approach reduces to solving a set of method's characteristics.
linear equations (Eq. 9) with an overall behavior similar to Newton's We would like to close by thanking all of our discussors for their
method. Because constraints are not forced to become feasible in one interest in explaining, comparing and improving different aspects of the
shot, but only gradually, as the minimization process takes place, the paper. This effort has no doubt improved the paper.
conditions described by Mr. Wollenberg have not occurred. It should be
noted that some methods do force violated constraints into the feasible
region by assigning them the limiting values together with some changes
of the variable set [ 16]. However, this is not the case in our method. In REFERENCES
fact, if constraints are such that it is impossible to satisfy all of them,
the method would stop at the best possible point.
The last point made by Mr. Wollenberg refers to the cases pre- [15] J. F. Dopazo, 0. A. Klitin and A. M. Sasson, "An improved initiali-
sented in the paper. The cases we picked to include in the paper just zation technique for the Hessian optimal load flow solution," sub-
happened to not have limiting inequalities at the optimum. However, mitted for the IEEE Summer Power Meeting, 1972.
violated inequalities did occur at the starting point (load flow solution) [16] J. Peschon, D. W. Bree and L. P. Hajdu, "Optimal solutions involv-
and these were brought in. In any case, as the method treats equality ing system security," Pager VI-B, 71 C26-PWR, 1971 PICA Con-
and inequality constraints in a similar fashion, and equality constraints ference Proceedings, pp. 210-218, 1971.

41

Authorized licensed use limited to: UNIVERSIDADE FEDERAL DE SAO CARLOS. Downloaded on April 29,2021 at 17:37:20 UTC from IEEE Xplore. Restrictions apply.