Identification of Load Characteristics of Power System Using GMDH

Copyri ght © I FA( : Itl t' lIl irica l iO I! ;l lId Sy... , t' 1l1 Estimation I m-!:~). York . l J K. 1~ I :-'G

Par;ll1H.' Ic..T

IDENTIFICA TION OF LOAD CHARACTERISTICS OF POWER SYSTEM USING GMDH S. Sagara*,

J.

Murata* and K. Kishimoto**

* 1)('/)(/1'/111 1'11/ of I': /i'l"/I'im/ FlIgilll'l'I'illg . !\.vll.l!ill U ll il' l'Isi/.\'. Flllwoka HI 2. j a/Jll II ** /? I'.II·a/'l'!i /) 1'/J11 I'/IIII'II/ . K \'II.I!i1l f;/I'(/ I'il l 'owl'/' ( ;0 .. l ilt .. Fllklloka HI 5, j a/){{ 11

Abstract. It is essential to power syste~ stability studies to represent load characteristics exactly. The proposed identification algorithms using GMDH produce , from field test data grouped into two sets , load representations iteratively in Taylor expansion form : at each iteration or layer , a term representing characteristics common to both data sets is selected and combined with the current representation to form a bette r one. The representation structure is of primary concern ; thus , the procedure, toge t her with the function form , is effective in preventing the load representations being complicated and in relating them to conventional l y used ones . Examples demonstrate the validity of the algorithms. Keywords. Identlfication; load characteristics .

modeling ;

GMDH;

INTROllUCTION In stability studies of electrlc power systems , represe n tations of load characteristics are necessar y . In early years of stability studies on digital computers , power system loads were rep r esented as constant impedances. However , acutual p ower system loads are composed of many devices which have various characteristics; as a result, system loads have rather complicated characteris tics different from those of constan t impedances . In the meanwhile , it has been po i nted out that load characteristics have much influence on the results of stability studies , see e.g . IEEE Committee Report (1973). Therefore a number of field tests (e . g. IEEE Committee Report , 1973 ; Nguyen , Ro b ichaud and Podesto, 1978) have been performed to get the data about load characteristics, and l oad representations have been constructed using the data . However, the presence of undesirable variation s i n the data mainly caused by consumers sometimes cause two kinds of probl ems : one is t h at parame t ers in load representations c h ange abruptly when the data change only sligh tly; the other is the contrary problem that the characteristics of different data can not be distiguished on their represen t ations. In this paper, algorithms for identification of load characteristics are proposed . The algorithms , using GMDH ( Group Method of Data Handling ) proposed by Ivakhnenko (1968) , suppress the problems mentioned above. The main feature is that function form of load representation is not fixed but adjusted to the nature of the data in the manner of h euristic self- organization , the idea descr i bed in the paper of Ivakhneko (1970). Thus they can extract effectively characteristics common to arbitrarily c h osen two sets of data . In other words, they can eliminate t he effects of compo nents which are not common to both data sets such as undesirable variations caused by consumers . Taylor expansion form is adopted as the framework , or the general form , of function forms of models; which is effective in preventing the models from being , as in ordinary GMDH's , too

power system ;

stability;

complicated a nd physically meaningless ones . Moreover the correspondences of the resultant l oad representations to those derived us i ng the conventional method i s clear. Examp l es that the algorithms are appl i ed t o some actua l f i e l d test data con f irm al l of a b ove mentioned merits .

INVESTIGATION OF THE PROBLEMS When one wishes to give a load representation, as its desirable features , high accuracy and/or ability to represent distinct character istics of dif f erent data, rather complicated representations are convenient . On t h e other hand, generality and small computation load, which a r e desirabl e features as wel l, prefer relatively s i mple ones. There must be a trade-off. In conven tional method , t h e f orm of loa d representa t i ons i s fixe d to s u ch as use d by Nguyen , Robichaud and Podesto ( 1 978): P Po

( ~ (P

(la)

(:or ( ~t

(lb)

(V: ( P q

~ °0 p

~ °0

(2a)

( :o(P

Po

=

(:or

(2b)

where P , Q, V and F are real power , reactive power, voltage and frequency respectively; Po' 00 ' and

are

their

standard

values;

and

np , nq , fp , fq are constants. The fixed function form may be too complicated for some data , while for the other data it may be too simple ; actually this

S. Sag-a ra . .J. Mur-a la and K. KishilllOIO

304

is the reason of the problems attaching to the conventional method. In the following sections, we develop the algorithms in which the function form or the complexity of load representations is not fixed but adjusted to the nature of the data.

GMDH GMDH is a cybernetic method for modeling inputoutput r e l ations proposed by A.G . Ivakhnenko in the late 1960's; its basic algorithm has the following features ( see Ivakhnenk o , 1968; Ivakhnenko, 1970 ; Ivakh nenk o , 1971 ; I keda , 1 979 ): (i) I t has multilayered structure; (ii) In each layer, variables are combined each other to form " partial models " as candidates for "perfect model" of the output variable; (iii) From the resultant partial models , some useful ones are se l ected usi ng a criterion and pass t h r ough to t h e next layer to be combined each oth er ; (iv) To assure the uniqueness of the perfect model and its continuity on the data changes, some regularizat i on is performed using the information conveyed in the data.

form . For each data set , the specific form is determined from t h e information contained in the data; the right hand side of each equat i on would be a finite series of some finite degree, and among the terms of degrees less than t h e maximum degree some might be omitted . Eqns. (3a) , (3b) , and (4b) can b e interpreted as Taylor ( 4 a) expansion form representations of static characteristics of load . This interpretation mak es it easy to take some physical meani ngs out of the representations . I f eqns. (la) and (lb) or (2a) and (2b) have enough accuracy , and if all difference values indicated by " 6 " in (3a) and (3b) or (4a) a n d (4b) are small, c ll and d l l in (3a) and (3b) wi ll b e nearly equal to np and nq in (la) and (lb), and cl and c2 in (4a) and (4b) will corresspond to np and nq in (2a) and (2b) respectively. Therefore eqns . (3a) , (3b) , (4a) and (4b) have correspondences to conventinally used representations such as eqns. (la) , (lb) , (2a) and (2b ) . Moreover they have abi 1 i ties fo r represent i ng characteristics more precisely.

Owing to these fea t ures , once criteria and algorithm structure wh ich are suitabl e for mode l ing purpose are given, an appropriate model structure wil l b e determined heuristically from the data , and a final model with the accuracy and the comp l exity fit for the purp ose will be der i ved . Th i s is just suitab le for sup p r essing the prob le ms in the mann er described in t h e preceding section. However t h ere is o n e th i ng specia l attention should b e payed to: in ordinary GMDH's the resul tant models are prone to be too complicated to tak e any physical meanings out of them.

- B 0,

k

_ A

~(k) Po

IDENTIFICATION ALGORITHMS -

B

IlP(k) Po

General Form of Load Representations

:§E(kl Po

0,

o

0

(e~Ar

-t rpA a k) p(k-l) + btk) 0 : estimated from Data A

a~ k) w PO(k-l) B

.

~ (6~r

+ b k) Vo

estimated from Data B

:

To preve n t load rep r esentations from being phys i ca lly meaningless ones and to avo i d tro u bles i n numer i cal calcu l ations for the data wh ich have only small variations , load representations are represented in the followi ng for m:

IlP(O) Po

w

+b(k) (e~r

a(k) p( k-l) 0 estimated fr o m the whole data

<

e > _Y _e:..s=--_ _ _ _--l

IlV) + c12 (IlF IlV ) 2 FO) + c2l ( Vo c l l ( Vo

+ c22

(e~ )( ~: )

+

C23 (~: )

2 +...

(3a)

e : (3b)

Po

cl

(e~ )

+

110 00

dl(e~ )

+

IlP

C2 (e~ )2 d2 (e~ )

+

2 +

C3 (e~ )3 d3 (e~)

Threshold Load repre s entation at t.h p k - th layer

: Final load representation +

3 +

...

(4a)

...

(4b)

where 6P=P-P ' 6V=V- V and 6F=F-F ; and 60 =0- 0 , O O O 0 ci , cij , d i and d i j (i ,j =l,2,3 ... ) are constants . Th e equations describe the framework , or ge neral

Values with superscript A, B or with no superscripts are c alculated or estimated from Data A, Data B or the whole data, respectively.

Fig. 1. Flow-chart of Algorithm V.

Identification of Power System Using GM OH

305

Method for Combination'

Meth od for Model Selection

Eqns . (3a), (3b), (4a) and (4b) are rather simple compared to other model structrues encountered in ordinary GMDH's. The simpleness prefers the method for constructing the mo del in which one variable in one layer is combined with the temporary partial model, like the method termed "successive determination" proposed by Ivakhnenko , Todua and Chukin (1972). Here, the term of the lowest degree is taken into the model first, then the terms of higher degrees are taken one by one. The presumption that a term of higher degree has less contribution and more sensitivity to the undesirable variations than a lower degree term and reduces comput ation speed in stab ility studies motivates such combination method.

Model selection criteria based on numerical accuracy of models are, though they are commonly used in ordinary GMDH's, likely to produce models with no generality but excessive complexity. A model selection criterion based on model structure should be developed . Such a criterion should be related to the components of models rather than the models themselves . As a preparation, the whole data are grouped into two sets . Then, t he terms which represents common characteristics of b oth data sets are selected, and the representations are composed of them. More specifically, coeffic i ent value of each term is estimated from each of two data sets, and the ratio of derived two estimates is calculated. If two estimates are nearly equal each other, i.e. the ratio of them lies around 1, the term can be regarded to represent the common c haracteriti cs . Terms with ratios within the predetermined threshold range are selected .

Algorithms

A -r.J:(O)

Po

~B

= 0, 1,

atk)

~p(O)

0,

j

~~k-1)

+

o

0,

Po

k

6tk)(e~)i-j ( ~:A) j

: estimated from Data A

( ~vB ) i-j ( F~FOB ) j

B B W B a\k) P (k- 1) + 6\k) Vo O : estimated from Data B a(k )

~:(k-l) + 6(k) (e~ rj

U:t

Flow- chart of an identification algorithm is shown in Fig. 1. We call it Algorithm V, because only voltage dependence of load is consi dered : load representations are produced to be of the f orm of eqns. (4a) and ( 4b ). The form is reasonable when the frequency dependence of load does not clear ly appears in the data. The whole data are grouped into two sets, A and B. The ratio of two estimates A B of coefficient value 6 (k)/6 (k) is compared to a threshold value 6 . If the ratio is within the A B threshold range ( 1/6 < 6 (k)/6 (k) < 6 ), the a lgorithm continues, otherwise it terminates . When frequency terms are considered, there will be produced several terms of the same degree. There

estimated from the whole data

150 3

;lE

0..

145 140

-1 0 No

... <1l

-1 5

>

;lE

0-

-20 - 25

No

1 20

e :

Threshold

> >

Final load representation

110 105

W(k) : Load representation at the k-th layer

Po

Value s with superscript A, B or with no superscripts are calculated or esti mated fr om Data A, Data B or the respectiv ely ,

115

whole data,

Fig. 2. Flow-chart of Algorithm F.

60 . 2 60 . 0 "59 . 8 N

:J:

0

50

100

150

200

ti me (s)

Fig. 3. Example of field test data.

250

S. Sagara,.J. Murata and K. Kishimoto

30!)

are no a priori knowledge about the superiority or inferiority among these terms. If the ratio is out of the threshold range, the next term will be combined and tested, until a final terminating condition meets. Figure 2 shows the flow-chart of the algorithm which produces representations of the form of eqns. (3a) and (3b), named Algorithm

Example of Application - 1 In data plots in Fig. 3, there are several steps which correspond to the tap positions. The data sets A and B are constructed in such a way that each step of part I in Fig. 3 is devided into two devisions equally.

F. Table 1 shows the coefficient values of the load representations derived from the data sets. For each load representation, residuals ( root mean square values) are calculated using part I,ll and III respectively, and plotted in Fig. 4. To make a comparison, the results of the conventional method which constructs eqns. (2a) and (2b) are also shown.

EXAMPLES OF APLLICATION TO FIELD TEST DATA

Data Data are derived from the field test in which the secondary voltage of the transformer in a substation is varied by changing the tap. Secondary real powel' P, reactive power Q, vol tage V and frequency F are recorded as averaged values over 6 cycles, or 0.1 seconds. Example of data is depicted in Fig. 3. Reactive power is compensated by synchronous condensors and static compensators equipped in the network, the compensated amount is estimated and added to the measured value.

TABLE 1

From Table 1 it can be seen that the coefficients corresponding to cl and d1 in (4a) and (4b) are nearly equal to those corresponding to np and nq in (2a) and (2b) respectively. The optimal value of threshold lies around 1.6 because using such threshold the resultant representations are relatively simple though they have enough accuracy. From Fig. 4, the accuracy in part I which is used for identification is improved as compared to the conventional method. The representations derived

Coefficient Values of the Load Representations from Part I (a) Real Power

algorithm

e

6V/V

Algorithm V

1.3 1.6 1.9 2.2

1.00 1.00 1.00 1.00

-1.05 -1.05 -1.05 -1.05

e

(~~)

U:l .. ·

1.3

1.02

Algorithm F

O

(6V/V )2 O

1.06

5.40

1.9

1.06

5.39

The conventional method

1.06 -

_4.0lx10 1.19x10 1.19x10

5.39

(~:(

... 4

1.15x10

...

(~~)U:r

6

3

2

1.48x10

2

1. 48x10

11 11

1.00

1

P

·. (~:t -1. 02x10

1.6

2.2

(~~)( ~:)

= 9.99x10- (V:l

Po

(b) Reactive Power algorithm

Algorithm V

Algorithm F

e

6V/V

1.3

6.14

O

(6V/V )2 (6V/V )3 (6V/V )4 (6V/V0)5 (6V/V )6 (6V/V )7 (6V/V0)8 (6V/V0)9 O o O O O 3.87x10

2

4.06x10

1.6

6.14

3.87x10

4.06x10

1.9

6.04

3.81x10

3.99x10

2.2

6.03

3.81x10

e

(~~)

.. .

1.3

6.43

4.60x10

1.6

6.33

3.99x10

2

3.99x10

(~~)

2

.. .

2 2

4.03x10 4.03x10

(~~t

6.06

3.82x10

4.01x10

2.2

6.06

3.82x10

4.01x10

!l Q

O

6.26 =

1.01(v:l

3

1. 30x10 1. 30x10

...

5 5

1.72x10 1.72x10

(~~f 9.88x10

1.9

The conventional method

3

2 2

4.04x10 4.04x10

6 6

.. .

4.12x10 4.12x10

7 7

(~~t

. ..

2 3 3

6.72x10

8

1.19x10

(~~ (( ~:) -2. 55x10

1. 31 x10 1. 31x10

5 5

7

1O

307

Identification of Power System Using CMDH by Algorithm V, though t hey are of s i mple forms, have no l ess accuracy than those derived using Algorithm F . Th e reason is that , as indi c ated in Fig. 3 , the freqency variation is not large.

proposed a l gorithm s are able t o represent the di s tinct charac teristics of each data set more c learly than the conventi ona l method.

The residuals f or part I and III are smaller than that for part 11 , which indicates , considering Fig . 3, t hat t he responses of load are different depending on whether vo ltage is being raised or being reduced. To make sure o f this, the residual s of the representations derived fro m part 11 us ing Algorithm V with the threshold 1 . 6 and using the conventiona l method are shown in Fig. 5. In this case , residuals for part 11 are s mal l while those for pa rt I and I I I are large. Besides , these differences are c l earer in the results of the proposed a l gor i thm than in the results of the conventional method . These confirm that the

Example o f Applicati on - 2

o

Data part Data part Il Data part III

Algorithm V Algorithm F The conventional method

•

o

0 . 75 ~

0.60

v

::l .-<

'"

> 0.45 (fJ

~ .-<

--G--

....'"

.-0-

~ 0 . 15

'"

Part I and III in Fig. 3 a re us ed as data A and B respectively. Coefficien t values of the l oad representa ti ons a r e s hown in Table 2 , and the residua l s a r e p l otted in Fig . 6 . Number s of terms in the representations are smaller t han in Table 1. The residuals for part I and III are closer each other than in Fig. 4; particularly, in the case of reactive power each val ue of both residuals and the difference between them ca l culated from the result of Algorithm V are sma l ler than those derived fr om the c onvent i onal method . Figure 7 shows the residuals of the load representation f r om part I, the data when voltage is being reduced , and part 11 , t he data on voltage being raised. The residuals for part I and 11 a r e very c l ose each other, and the residual for pa rt III is also small. These tel l us t hat these l oad represen tations re present common char acteristics of the data on voltage being reduced and those on voltage bei ng raised . The resu l ts described above shuw that the algorithms are more effecti ve than the conventional meth od in constructing, suppressing t he effect of undesirable variations, the l oad representations which are va l id over wider time range , and i n constructing load r epresentations wh ich can be used whether voltage is being raised or being reduced .

0.30

::l "tl

When preparing two data sets, there is an alternative way to that used in the pre c ed1ng subse c ti on: two of arbitrary data se t s a re chosen to be data A and B.

-0-

0.00 1.9

1.6

1.3

1. 0

Algorithms for identification of load characteristics of power system using G~1DH are proposed. Us ing the algorithms, l oad representat i ons with the f o llowing feat ures are obtained: it is easier than using the ordinary GMDH's t o interpret the

(a) Real power

Data pa rt Data part II Data part II I

CONCLUSIONS

2.2

Thresh o ld

Algorithm V Algorithm F 0 The conventional method

0

•

Data par t Data part II Data part III

0.25 0.20

..,"::l

'"

> Ul

.,,' ~

0 . 15

..,

'" ....

0.10

:..::s=

"tl

'""

- __

~

_____ O

-----e

;..

~:

e-

0.05

0 . 60 v

::l .-<

--D- -

::l

.,

.... .....

"~ • • • -O""

lE

'"

;If'

8

..

0

"......

•

:~

•

•

'"

1.6 )

_.0-

-<>0 .1 5 -0-

-0-

-<>-

_.0-

_.0-

-G-

-0-

(fJ

lE

'"

0. 30

0.1 0 - 0-

.-<

'"

::l "tl

....

0 . 00

;

_.0-

0 .45

>

0 Algor i thm V ( e 0 Th e conventional me thod

0 . 15

-0-

0.05

Ul

1.0

1.3

1.6

1.9

2.2



a:

0.00

0.00

Thresh o ld (b) Re a c tive power

Fig. 4 . Relationship between threshold value and residual ( part I ).

(a) Real powe r

(b ) Reactive power

Fig. 5 . Residua l s f or l oad representations ob t ained us ing the proposed a l gori t hm and the conventional method ( par t Il ) .

308

S. Sagara, TABLE 2

J.

Murata and K. Kishimoto

Coefficient Values of the Load Representations from Part I and III

physical meanings of them; they are ready to be taken correspond~nces with those from the conventional method and they can represent more detail. With the algorithms, the problems which occurred in the conventional method are suppressed. Moreover they can extract the characteristics common to two of arbitrary data sets; it is suitable for various purpose, e.g. constructing load representation which is valid in both cases of voltage being reduced and being raised.

(a) Real Power algorithm

6V/V

6

Algorithm V

1.3

9,77x10

1.6

9.77x10

1.9

9.77x10

2.2

9.77.10

6

Algorithm F

(6V/V0)2

O -1

1. 69

-1

1.69

-1

1.69

-1

1.69

...

6V/V0

(6V/V0)5

1.3

1.01

-5.05x10

1.6

1.01

_5.05x10

1.9

1.01

-5.05x10

2.2

1.01

_5.05x10

4

Algorithm V

Algorithm F

4

P V)9.76X10 Po = 1.00 (Vo

6V/V

1.3 1.6 1.9 2.2

6.36 6.36 6.36 6.36

6

6V/V

1.3 1.6 1.9 2.2

6.37 6.37 6.37 6.37

Data part Data part II Da ta part III

-0-

!l

QO

o

o

=

O

....'0'rn"

0.10

-0-

-0.-<

.. ( ) -

:.S:

~

::€r-

0.15

0.05

0.00

0.00

..().

..( ) -

(b) Reactive power

O

3.45x10 3.45x10 3.45x10 3.45.10

Fig. 7. Residulas for load representations obtained using the proposed algorithm and the conventional method ( part I and I I ).

O

REFERENCES

1. 01

(r .'!.-

26

Vo

Algor i thm V ( 6 = 1. 6 ) The conventional method

0.15

0.10

0.15

0.05

(a) Real power

0.30

~

(a) Real power

(6V/V )2

-0·

L-_ _ _ __

'"

>

Q)

0.30

0.00

.-<

Ul :E 0::

0.15

..

:l

0::

0.60

0.45

Algorithm V ( 6 = 1.6 ) The conventional mehtod

:l

6

The conventional method

0.45 Q)

4

(b) Reactive Power algorithm

0 0

4

-1 The conventional method

Data part I Data part II Data part III

0.00

~

-0-

::8= -<>-

•. 0-

L-_ _ _ __

(b) Reactive power

Fig. 6. Residuals for load representations obtained using the proposed algorithm and the conventional method ( part I and III ).

IEEE Commitee Report (1973). System load dynamics - simulation effects and determination of load constants. IEEE Trans. Power Apparatus and Systems, PAS-92, 600-609. Ikeda, S. (1979). Group method of data handling towards a modeling of complex systems. System and Control (Japan), 23, 710717 ( in Japanese ). Ivakhnenko, A. G. (1968). The group method of data handling a rival of the method of stochastic approximation. Soviet Automatic Control, 13, 43-55. Ivakhnenko, ~ G. (1970). Heuristic self-organization

in problems of engineering

cyberne-

tics. Automatica, 6, 207-219. Ivakhnenko, A. G. (1971). Polynomial theory of complex systems. IEEE Trans. Systems, Man and Cybernetics, SMC-1, 364-378. Ivakhnenko, A. G., M. M. Todua and Vu. V. Chukin, (1972). GMDH algorithm with successive determination of polynomial trends using the most essential variables. Soviet Automatic Control, 2, 44-54. Nguyen, C. T., Y. Robichaud and B. Podesto (1978). Determination of power-system load characteristics using a digital data-aqusition system. CIGRE, 1978 Session, 31-03.