678 IEEE Transactions on Power Apparatus and Systems,VVol. PAS-97, no.

3, May/June 1978

STATE ESTIMATION IN POWER SYSTEMS : DETECTING BAD DATA

THROUGH THE SPARSE INVERSE MATRIX METHOD

F. BROUSSOLLE
ELECTRICITE DE FRANCE

ABSTRACT In fact, using this formula makes it possible to keep at one's

disposal, following error detection, the jacobians and covariance
In the real time control of power system, the state estimation matrices calculated before error detection, hence to avoid a double
function plays a leading role. To perform its task, the estimator must calculation of these matrices.
detect and eliminate the grossly erroneous measurements.
The report begins with a summary of the principles of the
This report shows that one may take advantage of the sparse method. The second part deals with the calculations to be performed
structure of matrices in statistical computations for error detection for test and error detection. The last part is devoted to the analysis
tests. The sparse inverse matrix method makes it possible to perform of the results and performances of a program which applies the
statistical tests over hundreds of measurements, in real time, with proposed method.
performances allowing implementation of a complete estimator
(with accurate error detection) for a network of 300 nodes on a
small industrial computer. 2. PRINCIPLES OF THE METHOD

1. INTRODUCTION 2.1 -The mathematical model

Over past years, power system state estimation has been the Given a power system with known impedance and topology
subject of abundant literature. The first publications (1, 2 and 3) characteristics, n the number of buses and m the number of
laid the bases for electrical and mathematical state estimation models telemeasurements, at a given time, the electrical state of the power
and they highlighted the advantages of the various approaches. network is defined by the following system of equations:
The algorithms proposed for evaluating the state of the power z = h (x) + v
system are naturally derived from the traditional load flow calcula-
tions. In (1), the authors use the Gauss-Seidel solution method; in E(vvT) = W
all other methods, direct matrix calculations are used along with
treatments well-adapted to large size sparse matrices. These methods where
which include ordered elimination and bi-factorization (4, 5 and 6)
have proved highly efficient in the field of network problems. z : measurement vector
x : state vector (phases and magnitudes of voltages)
The first real time state estimator have experienced a certain v : measurement error vector
number of application difficulties especially as regards the complex W: covariance matrix
problem of grossly erroneous telemeasurements. Investigations were h : non-linear function inferred from network equations.
then concentrated on this subject and several authors (7, 8) have
proposed error detection methods. They showed that statistical tests
gave satisfactory results; unfortunately until now, this test required Assumptions on the modelling of errors
lengthy matrix calculations and it seemed impossible to use it in real
time for large size networks (exceeding 100 nodes). This situation - the mean error is equal to zero
was highly regrettable since the bigger the power system, the more - errors are small
numerous the risks of error. - errors comply with the normal distribution laws
- errors are non-correlated (matrix W is diagonal).
This report suggests a method which enables the statistical test
to be performed much more rapidly. It consists of taking advantage The measurements which do not meet the above-mentioned
of the sparse structure of the matrices in implementation of the tests, assumption are considered as grossly erroneous and should be remov-
by making use of the sparse inverse matrix technique. This technique ed. These erroneous measurements are detected and identified as
has been proposed by (9) to increase efficiency of the short-circuit shown in paragraph 2.3.
calculations which resulted in noticeable savings of time and storage
space. It makes possible the detection of errors, in real time, even in
the case of large size power systems. 2.2 - Power system state estimation ( 12)
Moreover, it should be emphasized that it is possible to combine Several algorithms can be used to estimate the state of a power
the sparse inverse matrix technique with the SHERMAN-MORISSON system. The most common and simplest one is the weighted least
formula which also permits computer time savings, as shown in (11). squares algorithm which provides very satisfactory results. It
consists of calculating vector x for which the following expression
is minimum:
J(x) = vTW l v
x can be calculated as the solution of the following equation

F 77 579-6. A paper recomrended and approved by (a J) = HTW 1 (z-h(x)) = 0

the IEEE Paier System Engineering Committee of the
IEEE Pcaer Engineering Society for presentation at the where H is the jacobian of h (x)
IEEE PES SurTer Meeting, Mexico City, Mex., July 17-
'2, 1977. Manuscript submitted January 17, 1977; nade The method to be used for solving this expression is similar to
available for printing April 6, 1977. that applied to a classical load flow calculation the NEWTON
0018-9510/78/0500-0678$00. 75 @)1978 IEEE
 679
method. It consists of transforming the non-linear system into a b. a search for the maximum value of ri/(F2 d' in order to locate
series of linear systems having the following form: this error. 1
xk+1 = xk + [HTw-1H]1 HTW-1 [Z-h(xk)] The first operation requires only the calculation of J(x) which
is both simple and rapid. The second operation requires the calcula-
where the series xk converges to the optimal solution k, if it exists. tion of the variances of estimation errors, which are the diagonal
elements of matrix E r:
The matrix HTW-V H is a symmetrical sparse matrix. Instead of
calculating its explicit inverse, it is more convenient to put it in a 1r = W-H E HT
factorized form, as follows
with
HTW -1H = LDLT
where: i =
[HTW- H]
Sections 2.2 indicates that for obtaining an optimum state
L Lower sparse triangular matrix estimation, it was not necessary to calculate explicitly the inverse
D Diagonal matrix Ex of HT W-1H. To take advantage of the sparse and symmetrical
LT Transposed fonn of L characteristics of this matrix, it is more efficient to put the
expression in factorized form, as follows:
Using this factorized form, at each iteration, xk+ 1_xk can be
calculated by solving the system HTW -1H = LDLT
LDLT (xk+ l_xk) = HTW1 (z-h(xk)) Unfortunately, to evaluate the variances, it may seem necessary
which results in significant time saving since direct inversion of to calculate the explicit inverse Ex of HTW1'H. However, at close
matrices can be avoided.
investigation, the explicit calculation of the inverse is not necessary.
The only elements which are indispensable are those used for calculat-
2.3 - Statistical tests and error detection ing a diagonal element of z r which can be expressed by the following
formula:
Let Sx = x-x be the error made on the state estimation of the
power system and the residue r = z-2 be the deviation between the
measurements and their estimated values; it can be shown that (E)dr2 - W
R
- j Lk Hi(HTW-1H)-1
R
kj HTQk
(8, 12):
E(Sx) = E(r) = 0 The elements (HTW H)1 k required are those corresponding
J
E(S x.S xT) = x = (HTW-'H)-l to products Hi HTQ # 0 which only involve elements other than
2 k
and that zero of matrix HTW-1H before inversion. Thus, the following rule
can be given:
E(r.rT) = r = W-HE xHT
,

"The elements of matrix Ex to be used for calculating variances are

Based on the error assumptions made in 2.1 concerning measure-
ments, the residues r are ruled by a zero mean normal distribution, those which correspond to the elements other than zero of the matrix
the variances of which are the diagonal elements of Er J(x) follows before inversion
a x2 law with an approximation of m-2n+ 1 degrees* of freedom. In addition, since the matrix HTW1'H is symmetrical, it is
possible to utilize the "sparse inverse matrix" method. This method
If certain measurements are grossly erroneous, the normality introduced a few years ago, (9) and (10), for the calculation of short-
assumptions no longer apply, r and J no longer follow the mentioned circuit currents is well adapted to the case of sparse matrices; it is
laws. The presence of a possible error is detected by a test on J and rapid and does not require substantial space in the memory.
the error is then located by a residual search on r.
If J(x) > x2 m-2n 1 p, one may deduce with a probability of 3.2 - Sparse inverse matrix method
error of p that at least one measurement is erroneous. The measure-
Given A, a symmetrical sparse matrix of which the inverse Z is
ment with the highest ri/(Er). is most likely the erroneous measure- to be calculated
ment and it is considered as a bad data. AZ= I,Iistheunitmatrix
Let us replace A by its factorized form
3. ERROR DETECTION CALCULATIONS
LDLTZ = I
3.1 -Introduction This relationship can be rewritten as follows:
An error is detected through two successive operations: D'L1'
LTZD
a. a test on J(R) for detecting the presence of an error. Put T = I - LT
Matrix I is a simple transformation of the matrix LT, the sign
* m is the number of elements of vector z, 2n-1 the number of elements having changed and the diagonal elements having been replaced by
of vector x zeroes. Z can then be expressed as follows:
680
Z= D'1L'1 + TZ The tests gave performance data of the new method compared
to the conventional method.
The organization of the program and its performances are
described in the sections below.

4.1 - Technical features of the program

Overall organization of the program

The following observations can be made:
The program consists of three main stages:
1. Matrix L is a triangular matrix with unit diagonal, its inverse has
similar properties. Thus, the diagonal of D-1 L-1 is equal to D-1 . a. Calculation and factorization of the jacobian
b. Calculation of the best estimation k
2. Matrix Z is symmetrical since A is symmetrical, therefore, to c. Calculation of the covariance matrix and detection of errors.
calculate Z, it is sufficient to calculate its upper part Zs3 i.e. the
ZJ elements such that j > i. Former variance calculation method
These two remarks show that matrix L-1 is not used in the The conventional program already contained a calculation
calculation of Zs' In fact, Zs complies with the following relation
sequence on measurement error variances. Matrix Ex = (HTW-lH)
zs/ D-1 + TZ
was calculated on a row-by-row basis. The contribution of each row
to the diagonal of Er was successively calculated. Matrix Ex was
therefore explicitly calculated but only the currently processed row
.EE was stored in the memory.

+ 1X~~~~~~~~~~~~ 4.2 - Performances of the program

It may be noted that owing to the peculiar form of T, the Size of the program
calculation of the elements of row k in Zs requires only the utiliza-
tion of the elements of the next rows. Z can be evaluated by The program can handle a system of the following maximum
successive calculations, starting with the last row. dimensions:
zJ=Ltk.zJ j>i, k>i 300 buses
i k i k 700 lines
2 000 telemeasurements (injections, flows and voltages).
zi= di+ E tkzi k >i for diagonal elements.
i i k i k Computation time
To calculate ZJ it is necessary to find first all the elements ZJ for The program was testes on a 250 node power system, with a
1 k measurement scheme consisting of 1200 telemeasurements. The
which k is such that tk is different from zero. This rule is also utilized computation times, recorded by an IBM 370-168 computer, yielded
the following results:
to create elements in L corresponding to the zero elements of A (in
addition to the elements in L which correspond to non zero elements a. Calculation and factorization of the jacobian 4s
of A). b. Calculation of the best estimate 4s
c. Calculation of variances and detection of errors
In order that the calculation may be completed, i.e. carried on - former method 15 s
until the first row of Z, it is sufficient to calculate the elements of Z - new method 5s
which correspond to non zero elements of L.
In addition, it is to be noted that the program is smaller when
In practice, on the programming level, the same core locations using the new method. Thus, thanks to the "sparse inverse matrix"
may be used for both L and Z, with the elements of Z gradually method, the program has become almost twice as rapid as the former
replacing the elements of L. one (13 seconds against 23 seconds), which represents a significant
time and memory saving.
When the calculation is completed, all the elements of Z
corresponding to elements of A other than zero have been calculated;
the matrix so obtained is referred to as the sparse inverse matrix of A. 5. CONCLUSION
This sparse inverse matrix contains only the basic elements of the
inverse matrix. In particular, it contains all the elements which occupy
the places of the elements other than zero in the initial matrix. More- In our opinion, the introduction of the "sparse inverse matrix"
over, the sparse inverse matrix makes possible an easy and direct method for the detection of errors may be regarded as a major
calculation of any element of the inverse whatsoever. progress. In addition, the application of the SHERMAN-MORRISSON
formula (11) makes it possible to proceed rapidly with the new
estimate following detection of one or several errors. Thanks to these
4. TEST RESULTS two techniques it has now become possible to obtain a state estima-
tion and error detection program sufficiently rapid and simple to be
The method described in this report has been extensively tested implemented in real time on an industrial computer to process systems
off-line on the 250 bus EDF high voltage system. of up to 300 buses.
 681
REFERENCES and model-decoupled. The algorithm-decoupled estimator converges to
the same solution point as the conventional least-squares estimator. The
model-decoupled estimator gives good approximations when the X/R
1. J.G. SIROUX and M. ADNET "Calcul en temps reel des modules ratios are high.
et des phases de tensions aux nceuds du reseau de transport A Besides the algorithm itself, we have applied detection and iclen-
partir des telemesures de puissances actives et r6actives". RGE tification in the decoupled mode. In detection, the index J is calculated
Vol 76 March-April, 1967, pp. 624-628. separately for the active and reactive power components of the pro-
blem, after which decoupled identification is performed. That is, if the
gross error detected is in the active power component, the identification
2. F.C. SCHWEPPE "Power system static state estimation" IEEE process is applied only among active power measurements.
Trans PAS-89 pp. 120-135, 1970. It seems to us that the decoupled estimation, detection and iden-
tification approach, using Z-sparse technique in the latter, offers con-
3. R.E. LARSON, W.F. TINNEY "State estimation in electric siderable improvement in the performance of the overall process, in
power systems" IEEE Trans PAS-89 pp. 345-363, 1970. terms of both computation and reliability, according to our tests on net-
works with high and low X/R ratios in the S.E. of Brazil.
4. M. CANAL, J. CARPENTIER "Ordinated elimination" PSCC We would like the authors' comments on the utilization of
LONDON 1965. generalized Sherman-Morrison formula, instead of retriangularization,
when the number of bad measurements is greater than one, especially
5. K. ZOLLENKOPF "Bi-factorization Basic computational about the maximum number of bad measurements that could be treated
-
with advantage.
algorithm and Programming techniques" Excerpt from "Large
Sparse Sets of Linear Equations" - ACADEMIC PRESS 1971
pp. 75-96.

6. W. TINNEY, W. POWELL, J. WALKER "Programming of sparsity

directed ordering schemes" PSCC CAMBRIDGE 1975, paper B. Magnus Manson (Power Division, ASEA Vasteras, Sweden): This is
3.1/1. an important contribution to the knowledge required for implementa-
tion of state estimation in energy control centers. Detecting and
7. H.M. MERILL et F.C. SCHWEPPE "Bad data suppression in eliminating gross measurement errors is vital in order to get a good
result from an estimation procedure using a weighted least squares
Power System State Estimation" IEEE PAS-90 pp. 2718-2725 criterion. The described method is elegant and seems to be efficient.
1971. However, it has some drawbacks. First. it mav be more difficult than
necessary to detect not too large deviations due to the smoothing effect
8. J.F. DOPAZO, O.A. KLITIN, A.M. SASSON "State Estimation of the estimation procedure, which makes several residuals deviate from
for power system detection of gross measurement errors" 8th zero due to one bad data.
PICA Conference MINNEAPOLIS June 1973. Secondly, it is rather timeconsuming to perform a new estimation
procedure once suspected bad data are located. Several reiterations may
9. TAKAHASHI, FAGAN, CHEN - "Formation of a sparse bus be required before a sufficient result is obtained.
impedance matrix" 8th PICA Conference MINNEAPOLIS A way to avoid these drawbacks is to search for gross measurement
June 1973. errors in the raw measurements, i.e. before state estimation. Because
any bad data detection procedure must make use of redundancy in the
measurement set - if there is no redundancy, there is no way to detect a
10. K. ZOLLENKOPF "Sparse nodal impedance matrix generated bad data - it is quite possible to detect most gross measurement errors
by the bifactorization method and applied to short circuit by simple checks on the raw data. We have found that checking, e.g.
studies". PSCC Cambridge 1975, Paper 3.1/3. nodal sums with far end measurements replacing unavailable near end
measurements, power flows in parallel branches, and power flows in
11. A. MERLIN - F. BROUSSOLLE "Fast method for bad data both ends of branches is a very efficient way to prevent most of the
identification in power system state estimation" Symposium gross measurement errors to reach the state estimator. Of course, it is
IFAC MELBOURNE 1976. beneficial to make an analysis of the residuals to detect any remaining
bad data. In this way, the desired result, a validated state estimate, may
12. A. LE ROY and P. VILLARD "Application of state estimation be attained more securely in a shorter time.
methods to evaluation of a telemeasurement configuration for The justification of the sparse inverse matrix method in the paper
energy power systems" excerpt from "computerized operation
seems to be a bit unfair to alternative methods. It is not necessary to
of power systems" published by Savu Crivat Savulescu - SAO
calculate the completely filled matrix X. The matrix H and the factors
Land D, which are all sparse, are already available and forward and
PAULO 1976. backward substitutions and matrix multiplications yield the desired
diagonal elements of L. This should be a considerably faster calcula-
tion than the complete inversion row by row followed by matrix
multiplications. Thus the figure 15 s may not be applicable for the best
alternative to the proposed sparse inverse matrix method.
Discussion Manuscript received August 9, 1977.

A. Garcia and A. Monticelli (UNICAMP, Campinas, Brazil): The

author is to be complimented for the paper presented.
There are two important points in this paper: the utilization of
sparse inverse matrix method and the application of the Sherman-
Morrison formula in the detection of bad data. F. Broussolle: The discussors comments and suggestions are most
The identification of measurements with gross errors is one of the welcome since they form guidelines as to how to improve state estima-
most important factors in static state estimation for electric power tion efficiency.
systems. However, the methods based on the calculation of the The question of Mm. Garcia and Monticelli deals with two very at-
covariance of measurements errors encounter a serious computational tractive techniques which are very useful to save time and memory in
limitation in the formation of the matrix of the covariance of the state estimation.
measurements. Without doubt, the use of Z-sparse calculation techni- Active and reactive decoupling, their first suggestion, gives to
que constitutes a significant advance in this area, by providing an effi- mean square state estimator the usual benefit that it gives to classical
cient method of calculating directly the diagonal elements of the load-flows: it reduces the iterative computations toward the solution R.
covariance matrix. Furthermore decoupling allows an additional improvement in the
Our experience in this area is based on the use of the Z-sparse tech- search for errors. It avoids searching the bad data in the reactive
nique as described in the paper, applied to fast decoupled state measurements when the error spoils only the estimated phase-angles
estimators. In this type of estimator, the active and reactive gain and reciprocaly searching in active measurements for an error in the
matrices remain constant throughout the least-squares iterative process. estimated modulus.
Two types of decoupled estimators were tested: algorithm-decoupled The second part of the question deal with the use of Sherman-
Manuscript received August 9, 1977. Manuscript received November 2, 1977.
682
Morrisson formula and Woodbury formula. This question has been the that it is always better to detect bad data before state estimation with a
subject of a paper (1) presented a few months before the present one. cheap algorithm on raw measurements than to detect them afterwards
The main idea of that previous paper is the following: with the heavy tools of statistical analysis. Simple checks on raw
Let us denote by HL. The matrix formed by the rows of H cor- measurements must be done whenever they are possible using local
responding to the measurements eliminated. L represents the set of in- redundancies in the set of measurements. Nevertheless, it must be notic-
dexes I of the measurements eliminated. HI is deduced from H by ed that simple checks as such as mentioned by Mr. Manson have a
removing the rows corresponding to these measurements. limited range of validity:
We can then write, using the Woodbury formula which generalizes
the Sherman-Morrisson formula: 1) checks assume linear modelisation of the network. They cannot be
applied to reactive and voltage measurements.
X7= IX + IXHL[SLL] HL Xx
2) very frequently the check will only tell you that there is one er-
with: roneous data among the active measurement of a substation and it
will not tell you which one is wrong.
SLL = WLL -HLHL
Likewise: With regard to the last question it must be precised that the pro-
I' I, + H' 7.H4:SLLV HLYxHIT gram with the figure of 15 s is not a rough one built only for the purpose
of giving a propitious reference to the new proposed method. It is based
The calculations are increasing with the number of measurements on careful sparse-oriented computations and is very close to the alter-
eliminated. One is led, furthermore, to inverse the matrix S1, whose native method of Mr. Manson.
dimension is equal to this number. However, up to more than 5
measurements eliminated, it is sure that Woodbury formula is a good
alternative to a new complete recalculation of E'r.
Numerous tests performed with active-reactive decoupling and REFERENCES
Woodbury formula on the EDF high voltage system gave excellent
results. 1. A. Merlin - F. Broussolle "Fast method for bad data identification
Addressing Mr. B. Magnus Manson's first suggestion about the in power system state estimation" - Symposium IFAC Melbourne
interest of bad data detection before state estimation, it must be said 1977.