Electrical power network modelling method

Information

  • Patent Grant
  • 6202041
  • Patent Number
    6,202,041
  • Date Filed
    Wednesday, April 29, 1998
    26 years ago
  • Date Issued
    Tuesday, March 13, 2001
    23 years ago
Abstract
A versatile modelling for any standard power system component is disclosed, by which the component modules can be easily plugged into the network module to form a small perturbation state space model of the entire system irrespective of the complexity of the system. The state equation is available as explicitly a function of every parameter and the input and output can be any variable, thus providing insight into the physical nature by simple matrix manipulation.
Description




BACKGROUND OF THE INVENTION




1. Field of the Invention




The invention relates to electrical power network modelling methods for small perturbation stability in a power transmission network.




2. Description of Prior Art




Power system stability has been and continues to be a major concern in the system or network operation. Large and small disturbance are in two main categories for stability analysis. In the past, the small disturbance performance was evaluated by the result of a transient stability program under a small perturbation. However, this simulation result will provide a limited insight because of the difficulty in taking measurements and in ensuring sufficient stability margins for swing oscillations. The small disturbance analysis is increasingly recognized because the spontaneous nature of swing oscillations can be analyzed based on a linearized system at the steady-state operating point. An eigenvalue analysis of the system described in this application can provide many insights which are difficult to be observed in transient plots.




Many methods have been proposed to represent networks, machines and associated control equipment such as the excitation system (EXC), governor system (GOV) and power system stabilizer (PSS) as well as new components of FACTs (flexible alternating current transmission) under small perturbation. In power system, the network and components are described by equations, and the control equipments are usually represented by blocks. In all other existing techniques, the control blocks are eventually transformed into equations in order to integrate with the network/component equations to form the state space equations. However, they have the following weakness:




(i) limited flexibility,




(ii) difficulty to interface with any user's new devices,




(iii) restricted input/output signal selection, e.g. ΔP


m


and ΔV


ref


for input signals,




(iv) infinite busbar assumption or restriction to a small hypothetical system,




(v) difficulty for computer program implementation,




(vi) limited exploitation of eigenvector analysis.




SUMMARY OF THE INVENTION




According to the invention there is provided an electrical system network and component modelling method for small perturbation stability in a power system, the method comprising converting the system network into elementary transfer blocks, converting the components into elementary transfer blocks, and plugging the components into the network to form a state space model of the entire system.




The invention is distinguished from the existing techniques that, instead of converting block to equations, all equations (network and component) are converted to blocks, such that blocks of control equipment (to any degree of complexity) can be easily is amalgamated by assigning simple node numbers. Because the network equations is also converted to blocks, any block modules of components, for instance machines, static var compensator (SVC), phase shifter (PS), high voltage direct current (HVDC) system and control and tieline, can be plugged into the network module.











BRIEF DESCRIPTION OF THE DRAWINGS




Embodiments of the invention will now be described by way of example with reference to the accompanying drawings in which:





FIG. 1

is an overall view of a Plug-in Modelling Technique (PMT) connection;





FIG. 2

shows elementary transfer blocks;





FIG. 3

is a network representation with one machine and one other component;





FIG. 4

is a third order machine module;





FIG. 5

is a fourth order machine module;





FIG. 6

is a fifth order machine module;





FIG. 7

is a sixth order machine module;





FIG. 8

is an SVC module;





FIG. 9

is a tieline module;





FIG. 10

is an example of a 4-block system;





FIG. 11

is an example of “non-elementary” transfer functions;





FIG. 12

is a formation of a ΔΩ


COI


I;





FIG. 13

is a tieline infinite busbar network module;





FIG. 14

is a configuration of an SVC;





FIG. 15

is an SVC module;





FIG. 16

is a general SVC voltage control;





FIG. 17

is configuration of a TCSC;





FIG. 18

is a TCSC module;





FIG. 19

is a configuration of a phase shifter;





FIG. 20

is a PS module;





FIG. 21

is a configuration of an AC/DC network at one terminal;





FIG. 22

is an equivalent T-circuit of a HVDC link;





FIG. 23

is a DC transmission line module;





FIG. 24

is an HVDC link module;





FIG. 25

is a tie line module;











DESCRIPTION OF PREFERRED EMBODIMENTS




Embodiments of the invention provide pragmatic and highly versatile methods, referred to as a Plug-in Modelling is Technique (PMT). Any standard power system components such as a static var compensator (SVC) a thyristor controllable series compensator (TCSC), a phase shifter (PS), a high voltage direct current (HVDC) system and control, a tieline, as well as a generator (of any degree of complexity) associated with control equipment can be modelled as modules and plugged into a network module (FIG.


1


). Because all these system components can be fully represented by two types of elementary transfer blocks and only five types of parameters as shown in

FIG. 2

(where “m” and “x” are the non-state and state variables respectively), the state space equation can be obtained easily by means of matrix manipulation irrespective of the complexity of the system. Stability of the system can be examined based on time/frequency response, eigenvalues, modal, and sensitivity analysis.




In the network system modelling, the concept of the PMT is derived from a previous Generalized Multimachine Representation (GMR). In GMR, however, loads are considered as constant impedances and the system components are restricted to machines. In the present application PMT, the network formulation is much enhanced by the inclusion of load characteristics and other system components. Moreover, options of different machine modelling are also provided.




For the stability studies, the load dynamics may be considered from:








S=P+jQ=P




o


(


V/V




o


)


a




+jQ




o


(


V/V




o


)


b


  (1)






Under small perturbation, the load characteristic can be represented by four fictitious admittances in (2).











[




Δ






I
R







Δ






I
J





]

=


[




Y
RR




Y
RJ






Y
JR




Y
JJ




]



[




Δ






V
R







Δ






V
J





]








where







Y
RR

=




(

a
-
2

)



PV
R
2



V
4


+



(

b
-
2

)



QV
R



V
J



V
4


+

P

V
2











Y
RJ

=




(

a
-
2

)



PV
R



V
J



V
4


+



(

b
-
2

)



QV
J
2



V
4


+

Q

V
2











Y
JR

=




(

a
-
2

)



PV
R



V
J



V
4


-



(

b
-
2

)



QV
R
2



V
4


-

Q

V
2











Y
JJ

=




(

a
-
2

)



PV
J
2



V
4


-



(

b
-
2

)



QV
R



V
J



V
4


+

P

V
2








(
2
)













and subscripts RJ stand for the network frame.




Busbar injections, voltages and network admittance matrix are related by I=YV i.e.










[




I
G






-

I
L





]

=


[




Y
GG




Y
GL






Y
LG




Y
LL




]



[









V
G












V
L





]






(
3
)













where subscripts G and L stand for generators and loads. Applying small perturbation and splitting the complex Y


NxN


matrix of (3) into real Y


2Nx2N


matrix (N=number of nodes), the load current I


L


is eliminated by equivalents given by (2), which is then included the block diagonal elements of Y, to obtain










[




Δ






I
G






0



]

=


[




Y
GG




Y
GL






Y
LG




Y
LL




]



[




Δ






V
G







Δ






V
L





]






(
4
)













Thus the load buses can be eliminated to











[

Δ






I
G


]

=


[

Y
GG


]





[

Δ






V
G


]









where






Y
GG



=


Y
GG

-




Y
GL



(

Y
LL


)



-
1




Y
LG








(
5
)













where Y


GG


′=Y


GG


−Y


GL


(Y


LL


′)


−1


Y


LG






Subsequently, only the generator buses are retained. However, if connected buses associated with other system Components are also considered, the network equation can be extended to










[




Δ






I
G








-
Δ







I
C





]

=


[




Y
GG




Y
GC






Y
CG




Y
CC




]



[




Δ






V
G







Δ






V
C





]






(
6
)













where ΔI


c


is the component current, ΔV


c


is the component busbar voltage and subscripts G and C stand for generator and other system components. Whilst generator terminals must be in series with a unique generator transformer, a component busbar can have multi-branch current flow. Therefore, the roles of the ΔV


c


and ΔI


c


in the network equation (6) must be swapped such that all components connected to the same busbar have a common input ΔV


c


. The modified network equation becomes










[




Δ






I

RJ
G








Δ






I

RJ
C






]

=


[




Y
GG




C
GC






C
CG




Z
CC




]



[




Δ






V

RJ
G








Δ






I

RJ
C






]






(
7
)













where Y


GG


′=Y


GC


Y


CC




−1


Y


CG


, C


GC


=Y


CC




−1


, C


CG


=−Y


CC




−1


Y


CG


and Z


CC


=Y


CC




−1


. Note that LHS (Left Hand Side) of (7), the “output” of network, is the “input” of components.




For instance, if the system consists one machine and one other shunt component, the full matrix equation will be










[










Δ






I
R1







Δ






I
J1











Δ






V
R2



_













Δ






V
J2





]

=


[



Y11


Y12


C13


C14




Y21


Y22


C23


C24




C31


C32


Z33


Z34




C41


C42


Z43


Z44



]



[










Δ






V
R1







Δ






V
J1











Δ






I
R2



_













Δ






I
J2





]






(
8
)













and these two components are connected by the network transfer blocks as shown in FIG.


3


.




In creating the representations of the network components, a machine is derived from Park's two-axis machine model. Machine behaviours can be described by different sets of equations of different order in the machine frame.




For instance, the 3rd-order machine equations in Park's dq-frame are








pE




q




′=[E




jd


−(


X




d




−X




d


′)


I




d




−E




q




′]/T




do


′  (9a)










pδ/Ω




0


=Ω−1  (9b)










pΩ=[P




m




−P




c




]/M


  (9c)










V




f




2




=V




d




2




+V




q




2


  (9d)










V




d




=−R




a




I




d




+X




q




I




q


  (9e)










V




q




=E




q




′−X




d




′I




d




−R




a




I




q


  (9f)










P




c




=V




q




I




q




+V




d




I




d


+(


I




d




2




+I




q




2


)


R




a


  (9g)






By applying small perturbation to the above equations:






Δ


E




q




′=[ΔE




jd


−(


X




d




−X




d


′)Δ


I




d


]/(1


+pT




do


′)  (10a)








Δδ=Ω


0




ΔΩ/p


  (10b)








ΔΩ=(Δ


P




m




−ΔP




c


)/


pM


  (10c)








Δ


V




r




=V




d




ΔV




d




/V




f




+V




q




ΔV




q




/V




r


  (10d)








Δ


V




d




=−R




a




ΔI




d




+X




q




ΔI




q


  (10e)








Δ


V




q




=ΔE




q




′−X




d




′ΔI




d




−R




o




ΔI




q


  (10f)








Δ


P




c


=(


V




q


+2


R




a




I




q





I




q




+I




q




ΔV




q


+(


V




d


+2


R




a




I




d





I




d




+I




d




ΔV




d


  (10g)






Equations (10) can be expressed by blocks or modules together with the block model of control equipment such EXC, GOV and PSS (FIG.


4


). Thus, a self-contained individual machine module is formed as shown in FIG.


4


. Other machine models can be similarly developed as shown in

FIGS. 5

,


6


and


7


.




For a machine and network interface, the generator will have its own dq frame which can be related to the network (RJ) frame by








I




d




=I




R


COS δ+


I




J


sin δ  (11a)










I




q




=−I




R


sin δ+


I




J


COS δ  (11b)










V




R




=V




d


COS δ−


V




q


sin δ  (11c)










V




J




=V




d


sin δ+


V




q


COS δ  (11d)






Under small perturbation, take the form:




 Δ


I




d


=COS δΔ


I




R


+sin δΔ


I




J




÷I




q


Δδ  (12a)






Δ


I




q


=−sin δΔ


I




R


+COS δΔ


I




J




−I




d


Δδ  (12b)








Δ


I′




R


=COS δΔ


I









d


−sin δΔ


V




q




−V




J


Δδ  (12c)








Δ


V




J


=sin δΔ


V




d


+COS δΔ


V




q




÷V




R


Δδ  (12d)






where δ is the angle between the dq frame and RJ frame.




Thus, the machine model can be “plugged” into the network by an adapter as shown in FIG.


3


.




Formation of other component modules (such as SVC, TCSC, PS, HVDC and tieline are provided in TFDEPC. All these modules, with inputs (ΔV


R


, ΔV


J


) and outputs (ΔI


R


, ΔI


J


), can therefore be “plugged” into the network module (FIG.


1


).




In network or system formation and analysis, using PMT the block structure is so general, computer algorithms (dealing only with two types of blocks and five parameters) can be easily programmed without a matrix format being concerned.




To run the program, a user has to




(a) prepare and perform a loadflow analysis




(b) input component data of the system (taking all the blocks, with their node numbers, within

FIG. 1

)




(c) input block data of the controller for the system components such as EXC, GOV and PSS (the “standard” models may also be retrieved from separate data files if preferred)




(d) start computation




The procedure of deriving the state equations (13) and (14), see below, or evaluating the sensitivity coefficients is unaffected by the complexity of the physical system is shown in TFDEPC 9.2 and TFDEPC 9.6. Because the entire system consists of only elementary blocks, existing software such as Matlab can be used as an analytical tool and advanced control techniques provided by the toolbox of Matlab can be applied directly and easily.








{dot over (X)}=AX+BR+E{dot over (R)}


  (13)










Y=CX+DR


  (14)






where the E{dot over (R)} term is often avoided in the usual representation but it does not create any analytical difficulty in the stability study.




Previously, the equipment representation (such as EXC) is limited to some specified models in order to save programming time for some “standard” packages. Due to the large variety of manufactured systems in practical use, simplification of this actual model to a specific model may provide a misleading result. It is very tedious to adjust the original EXC blocks (based on data easily obtained from the manufacturer's design) to fit them into some unrelated “standard” models. It would also be difficult to determine the optimum settings of specific control cards (e.g. the lead/lag time constants in the EXC feedback loop) since the “standard” model may contain no such blocks. Moreover, employing the “standard” model may lead to an improper PSS design, because the optimum setting is quite dependent on the EXC-PSS frequency response. The PMT of this patent application, however, facilitates the modelling of machines system components and their control equipment to any degree of complexity, including those IEEE models described by nonlinear functions or quadratic/polynomial transfer functions involving complex poles and zeros, as they can be converted readily to the elementary blocks of

FIG. 2

according to TFDEPC 9.4.




Due to the improvement in solid state technology, a lot of new components are introduced such as FACTS. In order to include these kinds of components, users need to form the plant using his own program because it is very difficult (if not impossible) for a general user to integrate these new components into the existing modelling methods. Due to the modular structure of the PMT, any new system component modules can be “plugged” in the network model via the specific input/output variables described in FIG.


1


. So there is no restriction to investigate system stability by consideration of the interaction among any components and the solid-state development of the power system networks applications can be enhanced.




In controller designs, the output vectors are restricted to state variables, e.g. speed or their linear combinations. This is probably because of the difficulty of deriving (14) for some non-state variable outputs e.g. ΔP


e


. In the so-called PQR techniques, the output variables are predefined to describe the network equation, whether they are used in subsequent system analysis or not. Due to this specified and limited input/output signals, it is difficult to take any unusual input/output signal for some typical controller design. In the PMT arrangement, any locally available variables such as voltage, current, MW, MVAr or MVA etc is can be used as control signals.




In the choice of reference frame, the infinite busbar, first introduced to reduce a system to a single machine equivalent but subsequently used to study multimachine systems, does not exist in reality. Instead of relying on an infinite busbar, an alternative assumption is that “network frequency” corresponds always to that of one arbitrarily chosen machine. However, in PMT, any machine speed can be used as a reference by means of a simple “−1” block (i.e. the connection of ΔΩ


ref


in FIG.


4


). A more reasonable or practical reference frame is to use the center of inertia (COI) of the entire network (or within a specific area) as shown in TFDEPC 9.5. Moreover, this COI frequency can be also used as a damping signal as described above.




Although eigenvalue technique has long been recognised as an effective tool to analyse power system small perturbation stability, the potential of eigenvalue analytical tools (modal and sensitivity analyses) has so far not been fully realized in other stability studies. The main restriction has been associated with the way in which the system state equations are normally formulated. The advantages of PMT method are brought out by the relative ease of its integration with the following techniques for practical synthesis.




For modal analysis in the PMT method, mode shape for any variable associated with the eigenvalue λ=α+jω can be determined easily from CU, where U is the eigenvector and C is the coefficient matrix in (14). Modal analysis for non-state variables was not applied. This may be attributed to the difficulty in establishing C if the y's in Y of (14) have to contain arbitrary non-state variables, for example, ΔP


e


, ΔI


d


or ΔI


q


. This is equivalent to expressing ΔP


e


, say, in terms of all the state variables (x's) of the machines, including EXC, PSS or governor, as well as of other components. In the PMT method, however, the y's can be chosen as any system variables or their combination without any restriction For sensitivity for arbitrary block parameters in the PMT method, matrix A is a function of all transfer block parameters so that the final expressions of ∂λ/∂κ are only simple algebraic scalar operations (e.g. ∂λ/∂b=wz/(hT


a


)) where the scalars (e.g. w,z and h) can be derived very easily from U and V as shown in TFDEPC 6.6.1. (“T


a


” and “b” are known transfer block parameters.)




The sensitivity expression for a change in block parameter is simple and straight forward. However, if the change is on a system parameter, numerous block parameters will be affected. For instance, the change in a bus MW loading P


Li


will affect all nodal voltages and the operating conditions of generators, which results in the change of some machine block parameters. Change in a line reactance X


ij


will furthermore affect the block parameters associated with the elements of matrices Y


GC


′, C


GC


, C


CG


and Z


CC


in the network equation (7). Fortunately, these changes, although very numerous, only affect the m-block parameters. With properly reordering the m-blocks as:






M=[M


1


, M


2


/M


3


]  (16)






where




M


1


=parameters varying with nodal voltages,




M


2


=parameters being the elements of matrices Y


GG


′, C


GC


, C


CG


and Z


CC


,




M


3


=the rest constant parameters




it is still possible to obtain ∂A/∂κ from the partial derivatives of the m-block parameters. Matrices Y


GC


′, C


GC


, C


CG


and Z


CC


associated with M


2


are created from the blocked matrix manipulation of nodal admittance matrix Y


GG


, Y


GC


, Y


CG


and Y


CC


in (6). While the block parameters belonging to M


1


, such as V


d


, V


q


, I


d


, I


q


, V


t


, cosδ and sinδ etc., can be expressed by busbar voltages (V and δ), even though these buses may have been eliminated in (7). Finally it is possible to evaluate ∂λ/∂κ, for κ=P


j


, Q


j


(arbitrary bus injection), G


ij


or B


ij


(arbitrary circuit admittance), because the differentiation of V or δ with respect to nodal injections has been determined in the loadflow algorithm.




The complicated sensitivity expression for system parameters are summarized in TFDEPC 9.6.2 for constant admittance load and third order machine model only. These equations can be directly applied to a higher-order machine model, or easily extended to other system components and load representations. What is needed is only the reordering of m-blocks and the reconstruction of M


1


, M


2


and M


3


. The sensitivity technique ∂λ/∂κ for κ=V


R


, V


J


has been applied to the probabilistic analysis of load variation , the sensitivity ∂λ/∂κ for κ=Q


L


has been applied to identify the best SVC location. The second order sensitivities ∂


2


λ/(∂κ


i


∂κ


j


) for κ


i


=V


R


, V


J


and K


j


=(PSS gain) can be used for evaluating the probabilistic distribution of sensitivity coefficient.




The validation of the PMT algorithms can be divided into two stages. The first stage (in three steps) is to compare the results with the well-known Heffron Phillips (HP) model.




(a) Machine modelling is first checked by using the new single machine infinite bus model as established in TFDEPC 9.8.




(b) The multimachine network described by (5) is checked against HP using a 2-machine 2-node system in which one machine (reference) has very large inertia and the output voltages ΔV


d


and ΔV


q


are made zero.




(c) This network equation (5) is once more checked by a large N-machine system in which two identical machines are connected to an arbitrary bus via identical reactance X


e


. (e.g. generator transformer reactance)




The above will confirm the validity of the 3rd-order machine model together with the network equation (5) because all the eigenvalues obtained in HP analysis reappear in all the above cases irrespective the complexity of EXC/GOV representations, although some more eigenvalues are present in (b) and (c) because there is more than one machine. The next stage is to check the network equation (7) (with any component) with the validated network equation (5) (with machines only):




(d) by treating an arbitrary transmission line (eventually eliminated in (5)) as a line component in (7),




(e) by treating a capacitor as a SVC component without feedback control,




(f) with the existence of both multi-SVC and multi-line components.




The validity of the network equation (7) can also be confirmed because eigenvalues obtained in computer runs are identical to the original (i.e. without SVC/line component). Note that the shunt SVC component connected to one node (

FIG. 8

) and the series line component connected between two nodes (

FIG. 9

) reflect the different ways of modular system connection.




Transfer Function Descriptions of Electric Power Components (Referred in the Specification as “TFDEPC”)




9.1.2.1 Static Var Compensator (SVC)




A typical SVC configuration with thyristor controlled reactor & fixed capacitor connected to a busbar through a step-down transformer is shown in

FIG. 14

where X


r


is the transformer impedance, B


C


and B


L


are the susceptance of capacitor and inductor. The SVC equations can be written as










X
SVC

=


X
T

-

1


B
C

+

B
L








(
108
)







[




V
R






V
J




]

=


X
SVC



[




-

I
J







I
R




]






(
109
)













Substituting eqn.108 in eqn.109 and applying small perturbation, eqn. 110 can be obtained and the block diagram of SVC is created in FIG.


15


. Besides, the general SVC thyristor controller block diagram is included in

FIG. 16

where V


SVC


and I


SVC


are the magnitude of voltage and current at the SVC terminal [16].










[




Δ






I
R







Δ






I
J





]

=


1

X
SVC




{


[




Δ






V
J








-
Δ







V
R





]

-



1


(


B
C

+

B
L


)

2




[




I
R






I
J




]



Δ






B
L



}






(
110
)













9.1.2.3 Thyristor Controlled Series Compensator (TCSC)




A typical TCSC configuration with thyristor controlled reactor & fixed capacitor connected in series with the transmission line is shown in FIG.


17


and described by










X
TCSC

=

-

1


B
C

+

B
L








(
111
)







[





V
R1

-

V
R2








V
J1

-

V
J2





]

=


X
TCSC



[




-

I
J1







I
R1




]






(
112
)







[




I
R1






I
J1




]

=

-

[




I
R2






I
J2




]






(
113
)













By applying small perturbation to the above equations and rearranging, eqn.114 can be obtained and the block diagram for the TCSC can be represented by FIG.


18


.










[




Δ






I
R1







Δ






I
J1





]

=


1

X
TCSC




{


[





Δ






V
J1


-

Δ






V
J2










-
Δ







V
R1


+

Δ






V
R2






]

-



1


(


B
C

+

B
L


)

2




[




I
R1






I
J1




]



Δ






B
L



}






(
114
)













9.1.2.4 Phase Shifter (PS)




Consider a phase shifter connected between busbars 1 & 2 and represented by voltage ratio T=1∠ψ in series with transformer admittance (y) as shown in

FIG. 18

, eqn.115 would be obtained.










[




I
1






I
2




]

=


y


[



1



-
T






-

T
*




1



]




[




V
1






V
2




]






(
115
)













By applying small perturbation to eqn.115, then










[




Δ






I
1







Δ






I
2





]

=



y


[



1



-
T






-

T
*




1



]




[




Δ






V
1







Δ






V
2





]


+


y


[





-
j







V
2


T






j






V
1



T
*





]



Δ





ψ






(
116
)













Splitting the complex matrix to real variable matrix, eqn.117 and

FIG. 19

can be obtained.










[










Δ






I
R1







Δ






I
J1










Δ






I
R2










Δ






I
J2





]

=



[




G
11




-

B
11





G
12




-

B
11







B
11




G
11




B
12




G
12






G
21




-

B
21





G
22




-

B
22







B
21




G
21




B
22




G
22




]



[










Δ






V
R1







Δ






V
J1










Δ






V
R2










Δ






V
J2





]


+


[










K
R1






K
J1









K
R2









K
J2




]


Δψ






(
117
)













i.e.









I




RJ




]=[Y




PS




][ΔV




RJ




]+[K




PS


]Δψ  (118)






9.1.2.5 High Voltage Direct Current System




Configuration of AC/DC network at each terminal is shown in FIG.


20


.




Converter Equation (117)








V




DC




=V




do


COS α−


R




C




I




DC


  (119)










V




R




I




R




+V




J




I




J




=V




DC




I




Dc


  (120)










I




2




+I




J




2




=k




DC




I




DC




2


  (121)






where







k
DC

=


18

π
2




a
2



k
2



n
2



B
2












V


do


=3(2)αV


t


B/π




and




is the communication overlap factor,




n is the number of poles per terminal,




B is the number of bridges in series for each pole.




a is the off-nominal turns ration of the transformer




Applying small perturbation,










Δ






V
DC


=




3


2


π


aB





cos





α





Δ






V
t


+


V
do


Δ





cos





α

-


R
C


Δ






I
DC







(
122
)







[




Δ






I
R







Δ






I
J





]

=


R


[




Δ






I
DC







Δ






V
DC





]


+

N


[




Δ






V
R







Δ






V
J





]







(
123
)













where






R
=



1
M



[






k
DC



I
DC



V
J


-


nI
J



I
DC







-

nI
J




I
DC









-

k
DC




I
DC



V
R


+


nI
R



V
DC







nI
R



I
DC





]


=

[




R
11




R
12






R
21




R
22




]






N
=



1
M



[





I
R



I
J





I
J
2






-

I
R
2






-

I
R




I
J





]


=

[




N
11




N
12






N
21




N
22




]












and




M=I


R


V


J


−I


J


V


R






Current Regulator Equation










p





cos





α

=




K
C


T
C




(


I
DC

-

I
ref

-

U
C


)


-


1

T
C



cos





α






(
124
)













Under small perturbation,










Δ





cos





α

=



K
C


1
+

p






T
C






(


Δ






I
DC


-

Δ






I
ref


-

Δ






U
C



)






(
125
)













i.e. The firing angle (α) in eqns.119,122,124 and 125 will become the advance angle (β) and I


DC


will become negative for inverter side.




DCA Network




If a T-model is used to represent the DC line (FIG.


22


), then








V




DC1




=R




1




I




DC1




+L




1




pI




DC1




+V




C


  (126)










V




DC2




=R




2




I




DC2




+L




2




pI




DC2




+V




C


  (127)

















V
C

=


1
Cp



(


I
DC1

+

I
DC2


)






(
128
)













Hence










Δ






I
DC1


=



Δ






V
DC1


-

Δ






V
C





R
1

+


L
1


p







(
129
)







Δ






I
DC2


=



Δ






V
DC2


-

Δ






V
C





R
2

+


L
2


p







(
130
)







Δ






V
C


=


1
Cp



(


Δ






I
DC1


+

Δ






I
DC2



)






(
131
)













Therefore, the block diagram of a DC transmission line is formed as FIG.


22


and the HVDC link module can be obtained as shown in FIG.


23


.




9.1.2.6 Tieline




Tieline flow is one of the very effective damping signal. In order to get/monitor this signal, the tieline (or any circuit) has to be treated as one of the system components. Since the terminal voltages & currents of a tieline can be related by:










[




I
1






I
2




]

=


[




Y
11




Y
12






Y
21




Y
22




]



[




V
1






V
2




]






(
132
)







[




Δ






I
1







Δ






I
2





]

=


[




Y
11




Y
12






Y
21




Y
22




]



[




Δ






V
1







Δ






V
2





]






(
133
)













By splitting the complex matrix to real variable matrix,






[ΔI


RJ


]=[Y


TL


][ΔV


RJ


]  (134)






and the tieline model is shown in FIG.


9


. Thus any tieline signals at either side can be obtained by










Δ





P

=



V
R


Δ






I
R


+


I
R


Δ






V
R


+


V
J


Δ






I
J


+


I
J


Δ






V
J







(
135
)







Δ





Q

=



V
J


Δ






I
R


-


I
J


Δ






V
R


-


V
R


Δ






I
J


+


I
R


Δ






V
J







(
136
)







Δ





S

=



P
S


Δ





P

+


Q
S


Δ





Q






(
137
)







Δ





I

=




I
R

I


Δ






I

R







+



I
J

I


Δ






I
J







(
138
)







Δ





V

=




V
R

V


Δ






V

R







+



V
J

V


Δ






V
J







(
139
)













9.2 Formation of State Space Equation




In the PMT, the blocks are either of zero order, as in

FIG. 2



a


, or of first order as in

FIG. 2



b


. The “in” and “out” of each block are assigned with node numbers. An example of connection of a four-block system is shown in FIG.


12


and referred to in TFDEPC 9.3. Thus, referring to the respective node numbers, the connection L-matrix can be constructed in the following form:










[







X
i





Y








M
i




]

=


[




L
1




L
2




L
3






L
4




L
5




L
6






L
7





L
8





L
9





]



[



X




R




M



]






(17a)
(17b)
(17c)














where X


i


, X, M


i


and M are vectors collecting all the x


i


, x, m


i


and m variables of

FIG. 2

, and R and Y are the input and output vectors. As the equation of the nth zero order block is m


n


=k


n


m


i


, the matrix equation for all zero order blocks will be






M=KM


i


  (18)






where








K=<k




n


>  (19)






=diagonal matrix of k


n


, n=1,2, . . . .




Similarly, the equation of the first order block is








x




n




=x




i






n




(


b




n




+pT




b






n




)/(


a




n




+pT




a






n




)  (20)






Expanding (20), the equation for all the first order blocks is







X=−K




c




X+K




b




X




i




+K




f




{dot over (X)}




i


  (21)




where








K




a




=<k




a






n






>=<a




n




/T




a






n




>  (22a)










K




b




=<k




b






n






>=<b




n




/T




a






n




>  (22b)






and








K




f




=<k




r






n






>=<T




b






n






/T




a






n




>  (22c)






are diagonal matrices (n=1,2 . . . ).




The state space equation can be constructed from (17), (18) and (21) by eliminating M


i


, M and then X


i


from (17).




Therefore, substituting (17c) in (18),








M=KM




i




=KL




7




′X+KL




8




′R+KL




9




M=L




7




X+L




8




R+L




9




M








so








M=H


(


L




7




X+L




8




R


)  (23)






where








H


=(


I−L




9


)


−1








and






I=identity matrix






Substituting (23) in (17b),








Y=L




1




X+L




5




R+L




6




H


(


L




7




X+L




8




R


)






so








Y=CX+DR


  (14)






where








C=L




1




+L




6




HL




7








and








D=L




5




+L




6




HL




8








Substituting (23) in (17a),








X




i




=L




1




X+L




2




R+L




3




H


(


L




7




X+L




8




R


)






so








X




i




=FX+GR


  (24)






where








F=L




1




+L




3




HL




7


  (25)






and








G=L




2




+L




3




HL




8








Substituting (24) in (21),








{dot over (X)}=−K




a




X+K




b


(


FX+GR





K




i


(


F{dot over (X)}+G{dot over (R)}


)






so that








{dot over (X)}=AX+BR+E{dot over (R)}


  (13)






where








A=S


(


K




b




F−K




a


)  (26)








B=SK


b


G








E=SK


c


G






and








S


=(


I−K




c




F


)


−1


  (27)






9.3 Formation of the L-Matrix




With reference to the simple example of a four-block system shown in

FIG. 12

, and by comparing the corresponding node numbers of the LHS and RHS vectors (“1” if equal otherwise “0”), the sparse connection-matrix equation is




















(
12
)






(
13
)









(
13
)









(
11
)









(
14
)






[













x

i
1







x

i
2












y




_

_









m

i
1










m

i
2





]


=



[



















1










1






























1












































1

















1











1



























(
13
)




(
14
)









(
11
)





(
12
)




(
11
)





]



[













x
1






x
2











r




_

_









m
1









m
2




]
















(
13
)






(
14
)









(
11
)









(
12
)









(
11
)















Since the input data are a sequence of blocks (i.e. branch) rather than nodal information, this technique is much simpler and faster than the indirect method of using incidence matrices in this case.




9.4 Representation of Complicated Transfer Functions




Nonlinear function under small perturbations, or “quadratic” transfer functions, may be expressed as elementary blocks in

FIG. 2

, as illustrated for some typical cases shown in FIG.


13


.







Case a:










Z
=





X





Y








Δ





Z

=






Y





Δ





X

+

X





Δ





Y











Case b:










E
=





&LeftBracketingBar;

V
+

j





I





X


&RightBracketingBar;








E
2

=






V
2

+


I
2



X
2










Δ





E

=







(

V
/
E

)


Δ





V

+


(


X
2



I
/
E


)


Δ





I











Similarly










I
=





&LeftBracketingBar;


I
d

+

j






I
q



&RightBracketingBar;








Δ





I

=







(


I
d

/
I

)


Δ






I
d


+


(


I
q

/
I

)


Δ






I
q












Case c:










E
=





f


(
I
)









Δ





E

=





m





Δ





I










where





m

=

slope






(



E

/


I


)






at





the





operating






point
.





Case d:    A 'quadratic' transfer function  may be
represented by two elementary blocks.  The  two block
diagrams are identical so long as
















G
=





C
1








A
=






D
0




C
1

/

C
0









K
=







C
0

/

C
1
2


+


D
0

/

C
0


-


D
1

/

C
1









B
=






C
0

/

C
1















In the limiting case where C


o


=0, C


1


=0, or if the numerator polynomial is also quadratic, it is still possible to express the original second order block as a combination of elementary blocks by more elaborate circuit block manipulation. Being a combination of quadratic and first order, the high order polynomial can be also represented by elementary blocks.




9.5 COI Reference Frame




The position of the COI is defined as










δ
COI

=




i
=
1

p




k
i



δ
i







(
28
)













where k


i


=H


i


/H


T


and H


i


, H


T


are the individual and the sum of the inertia constant of the generator. δ


i


is the position of generator i with respect to the stationary frame.










i
.
e
.





Ω
COI


=




i
=
1

p




k
i



Ω
i







(
29
)













Applying the small perturbation to (29), (30) and the COI speed signal are obtained (FIG.


14


).










Δ






Ω
COI


=




i
=
1

p




k
i


Δ






Ω
i







(
30
)













9.6 Evaluation of Sensitivity Coefficients ∂λ/∂κ




9.6.1. Sensitivity for Arbitrary Block Parameters




From (26) and (27), system matrix A is a function of K


b


, K


a


and K


t


. Hence, for κ=k


b


, k


a


or k


t


in the nth block as defined in (22), ∂A/∂κ can be expressed as follows:












A




k
b



=

S





K
b





k
b




F





(31a)










A




k
a



=


-
S






K
a





k
a










or




(31b)












A




k
t



=








S




k
t





(



K
b


F

-

K
a


)








=






-
S






(

I
-


K
t


F


)





k
t





S


(



K
b


F

-

K
a


)









=






S





K
t





k
t





FS


(



K
b


F

-

K
a


)



=

S





K
t





k
t




FA









(31c)













where ∂K/∂k (∂K


b


/∂k


b


, ∂K


a


/∂k


a


and ∂K


c


/∂k


c


) will have only one non-zero element of unity value (the nth diagonal element). Therefore, if W


T


=(w


1


, w, . . . ) and Z


T


=(z


1


, z


2


. . . ) are two arbitrary vectors, the product of W


T


(∂K/∂k)Z=w


n


z


n


will be as simple as the product of two scalars. Because of this feature, ∂λ/∂κ can be evaluated very easily as shown below.




Define the scalar






h=V


T


U=(v


1


, v


2


, . . . ) (u


1


, u


2


, . . . )


T








and the vectors






W=V


T


S=(w


1


, w


2


, . . . )  (32a)








Z=FU=(z


1


, z


2


, . . . )


T


  (32b)






where S and F are already available in (27) and (25), respectively, when forming the state space equation. Substituting (31) in (15),












λ




k
b



=




V
T


S





K
b





k
b




FU



V
T


U


=




W
T






K
b





k
b




Z



V
T


U


=



w
n



z
n


h







(33a)










λ




k
a



=




-

V
T



S





K
a





k
a




U



V
T


U


=

-



w
n



u
n


h









or




(33b)










λ




k
t



=




V
T


S





K
t





k
t




FAU



V
T


U


=




V
T


S





K
t





k
t




FU





λ



V
T


U


=



w
n



z
n


λ

h












Hence





for





κ

=
b

,
a
,

T
b

,

or






T
a


,




λ

/


κ







becomes

,

using


(
22
)


,





(33c)









λ



b


=





λ




k
b








k
b




b



=





λ




k
b





1

T
a



=



w
n



z
n



hT
a








(34a)









λ



a


=



-

w
n




u
n



hT
a






(34b)










λ




T
b



=



w
n



z
n


λ


hT
a








or




(34c)












λ




T
a



=









λ




k
b








k
b





T
a




+




λ




k
a








k
a





T
a




+




λ




k
t








k
t





T
a











=






-

(





λ




k
b




b

+




λ




k
a




a

+




λ




k
t





T
b



)


/

T
a
2








=






-


w
n



(



z
n


b

-


u
n


a

+

λ






z
n



T
b



)



/

(

hT
a
2

)









(34d)













Consequently, the sensitivities of λ, with respect to the nth block, are obtained by picking up nth elements from the vectors U, W arid Z and then performing a simple algebraic operation with some other scalars according to (34), where W and Z can be derived from V and U in (32). Therefore the sensitivity with respect to all block parameters can be evaluated once theses vectors are formed.




9.6.2 Sensitivities for Arbitrary System Parameters




With diagonal matrix K of (19) collecting all the parameters of zero order blocks, parameters associated with system operating condition are all included in K. From (25), (26), and (27), ∂A/∂κ can be expressed as












A



κ


=





[

S


(



K
b


F

-

K
a


)


]




κ


=

S


(



K
t





F



κ



A

+


K
b





F



κ




)







(35a)









F



κ


=





(


L
1

+


L
3


H






L
7



)




κ


=


L
3


H




K



κ




(



L
9



H

+
I

)



L
7








(35b)













Investigating the elements of K associated with system operating condition, some of them are the simple functions of variables associated with generator buses, which are collected in M


1


of (16) or M


1(i)


of (36) corresponding to i-th machine.








M




1(i)




=[I




di




, I




qi




, V




di




, V




qi




, V




ti


, cos δ


i


, sin δ


i


]  (36)






The others describing network equation are collected in M


2


of (16). With ∂K/∂κ constructed from the partial derivatives of elements of M


1


and M


2


, ∂λ/∂κ for κ=(system parameter) therefore can be evaluated from (15) and (35).




9.6.2.1 Derivatives of M


1(i)


to j-th Nodal Voltages (V


Rj


+jV


Jj


)




Firstly, assuming that the nodal voltages and generator currents are defined in rectangular coordinates. E


Qi


can also be expressed referencing network frame as:












E
QRi

+

j






E
QJi



=


(


V
Ri

+

j






V
Ji



)

+


(


I
Ri

+

j






I
Ji



)



(


R
ai

+

j






X
qi



)




,






E
QRi

=


V
Ri

+


I
Ri



R
ai


-


I
Ji



X
qi




,






E
QJi

=


V
Ji

+


I
Ri



X
qi


+


I
Ji



R
ai








(
37
)








E
Qi
2

=


E
QRi
2

+

E
QJi
2



,






E
QRi

=


-

E
Qi







sin






δ
i



,






E
QJi

=


E
Qi






cos






δ
i







(
38
)













From the nodal voltage equation of overall network I=YV, coordinates transformation (11) and,










V
ri
2

=


V
Ri
2

+

V
Ji
2






(
39
)













the partial derivatives can be expressed as follows by using a constant C


o


with C


o


=1 for i=j and C


o


=o for i≠j. Subscript “(.)j” stands for “Rj” or “Jj”.













I
Ri





V
Rj



=





I
Ji





V
Jj



=

G
ij






(40a)










I
Ri





V
Jj



=


-




I
Ji





V
Rj




=

-

B
ij







(40b)










E
QRi





V
Rj



=





E
QJi





V
Jj



=


C
0

+


R
ai






I
Ri





V
Rj




-


X
qi






I
Ji





V
Rj










(41a)










E
QRi





V
Jj



=


-




E
QJi





V
Rj




=



R
ai






I
Ri





V
Jj




-


X
qi






I
Ji





V
Jj










(41b)











E
Qi





V


(
.
)


j




=



-




E
QRi





V


(
.
)


j






sin






δ
i


+





E
QJi





V


(
.
)


j





cos






δ
i















sin







δ
i





V


(
.
)


j




=



-
1


E
Qi




(





E
QRi





V


(
.
)


j




+





E
Qi





V


(
.
)


j





sin






δ
i



)



,





(41c)











cos







δ
i





V


(
.
)


j




=


1

E
Qi




(





E
QJi





V


(
.
)


j




-





E
Qi





V


(
.
)


j





cos






δ
i



)
















I
di





V


(
.
)


j




=










I
Ri





V


(
.
)


j





cos






δ
i


+





I
Ji





V


(
.
)


j





sin






δ
i


+















I
Ri




cos








δ
i




V


(
.
)


j





+


I
Ji




sin








δ
i




V


(
.
)


j






,









(
42
)














I
qi





V


(
.
)


j




=







-




I
Ri





V


(
.
)


j






sin






δ
i


+





I
Ji





V


(
.
)


j





cos






δ
i


-














I
Ri




sin








δ
i




V


(
.
)


j





+


I
Ji




cos








δ
i




V


(
.
)


j





















V
di





V
Rj



=



C
0


cos






δ
i


+


V
Ri




cos








δ
i




V
Rj




+


V
Ji




sin








δ
i




V
Rj






,









V
di





V
Jj



=



C
0


sin






δ
i


+


V
Ri




cos








δ
i




V
Jj




+


V
Ji




sin








δ
i




V
Jj






,









V
qi





V
Rj



=



-

C
0



sin






δ
i


-


V
Ri




sin








δ
i




V
Rj




+


V
Ji




cos








δ
i




V
Rj






,





(
43
)










V
qi





V
Jj



=



C
0


cos






δ
i


-


V
Ri




sin








δ
i




V
Jj




+


V
Ji




cos








δ
i




V
Jj









(
44
)











V
ti





V
Rj



=


C
0




V
Ri


V
ti




,





V
ti





V
Jj



=


C
0




V
Ji


V
ti








(
45
)













9.6.2.2 Derivatives of M


2


to Nodal Voltages (V


Rj


+jV


Jj


)




M


2


of (16) will be consisted of the elements of Y


GG


′ alone when the machine is considered only. When static loads are represented by equivalent shunt admittances and added to the corresponding diagonal elements of Y


LL


′ of (4), only the equivalent admittances are affected by nodal voltages. Thus,
















Y
GG






V


(
.
)


j




=










V


(
.
)


j






[


Y
GG

-




Y
GL



(

Y
LL


)



-
1




Y
LG



]








=






Y
GL



Y
LL

-
1







Y
LL






V


(
.
)


j






Y
LL

-
1




Y
LG









(
46
)













When the deviation of nodal voltages is regarded being independent with load powers, differentiation of (47a) yields (47b).











Δ






G
ii


+

j





Δ






B
ii



=


(


P

l
,
i


-

j






Q

l
,
i




)

/

V
ti
2






(47a)










(


Δ






G
ii


+








B
ii



)





V


(
.
)


i




=




-
2



V


(
.
)


i




V
ti
4




(


P

l
,
i


-

j






Q

l
,
i




)






(47b)













If the deviation of nodal voltages is risen from the variation of nodal powers, the equivalent admittances should be expressed by nodal injected currents as:











Δ






G
ii


+

j





Δ






B
ii



=

-



I
Ri

+

j






I
Ji





V
Ri

+

j






V
Ji









(
48
)













and the injected currents are determined by nodal voltage equation I=YV. Therefore,














Δ







G
ii





V


(
.
)


j




=



-
1


V
ti
2




[



C
0



(

2


V


(
.
)


j



Δ







G
ii

÷

I


(
.
)


j




)


+


V
Ri







I
Ri





V


(
.
)


j




÷

V
Ji







I
Ji





V


(
.
)


j






]



,









Δ







B
ii





V
Rj



=



-
1


V
ti
2




[





C
0



(

2


V
Ri


Δ







B
ii

÷

I
Ji



)


÷

V
Ri







I
Ji





V
Rj




-


V
Ji






I
Ri





V
Rj





]



,









Δ







B
ii





V
Jj



=



-
1


V
ti
2




[



C
0



(


2


V
Ji


Δ






B
ii


-

I
Ri


)


+


V
Ri






I
Ji





V
Jj




-


V
Ji






I
Ri





V
Jj





]







(
49
)













9.6.2.3 Derivatives to Nodal Voltages (V


j


, δ


vj


)




When nodal voltages in polar coordinates are defined as V


j


∠δ


vj


, the derivatives of any element M


i


in and M


1


or M


2


of (16) can be obtained from V


Rj


=V


j


cosδ


vj


and V


Jj


=V


j


cosδ


vJ


as:














M
i





V
j



=






M
i





V
Rj




cos






δ
vj


+





M
i





V
Jj




sin






δ
vj




,









M
i





δ
vj



=



-




M
i





V
Rj






V
j


sin






δ
vj


+





M
i





V
Jj





V
j


cos







δ
vj

.








(
50
)













9.6.2.4 Derivatives of K to Nodal Injected Power (P


j


+Q


j


)




Jacobian matrix J used in load flow calculation determines the linearized relationship between nodal powers and nodal voltages. ∂V


Ri


/∂P


j


, ∂V


Ri


/∂Q


j


, ∂V


Ji


/∂P


j


and ∂V


Ji


/∂Q


j


are the elements of J


−1


. Therefore, the derivatives of any element M


i


in and M


1


or M


2


of (16) to nodal powers will be expressed as:














M
i





P
j



=



i



(






M
i





V
Ri








V
Ri





P
j




+





M
i





V
Ji








V
Jj





P
j





)













M
i





Q
j



=



i



(






M
i





V
Ri








V
Ri





Q
j




+





M
i





V
Ji








V
Jj





Q
j





)







(
51
)













9.6.2.5 Derivatives of K to Line Admittance (y


mn


=g


mn


+jb


mn


)




Change of a line admittance y


mn


will affect admittance matrix Y and the initial system operating condition described by nodal voltages and injected currents. Derivative of K with respect to b


mn


has the same form as that to g


mn


. Only the latter is discussed.




The derivatives of Y


GG


′ to g


mn


will be obtained from ∂Y/∂g


mn


which is constructed as: “1” on the sits of (m,m) and (n,n), “−1” on the sits of (m,n) and (n,m), “0” on other sites The derivative of Y


GG


′ with respect to g


mn


is:













Y
GG






g
mn



=





Y
GG





g
mn



-





Y
GL





g
mn






(

Y
LL


)


-
1




Y
LG


+




Y
GL



(

Y
LL


)



-
1







Y
LL






g
mn






(

Y
LL


)


-
1




Y
LG


-




Y
GL



(

Y
LL


)



-
1







Y
LG





g
mn









(
52
)













However, due to I=YV, (40) will become (in vector form)












I




g
min



=


Y




V




g
mn




+




Y




g
mn




V






(
53
)













From (41) to (45), if ∂M


i


/∂V is regarded as a linear function of ∂I/∂V denoted as f(.), the derivatives of M


i


will be represented in vector form as













M
i





g
mn



=



f


(



I



V


)






V




g
mn




+

f


(




Y




g
mn




V

)







(
54
)













∂V/∂g


mn


can be solved from (55) which is obtained by differentiating the nodal power equation.
















V
i





g
mn





I
i
*


+



j



(



V
i






Y
ij
*





g
mn





V
j
*


+


V
i



Y
ij
*






V
j
*





g
mn





)



=
0

,





i
=
1

,
2
,
3
,





(
55
)













9.7 Single 3rd-Order Machine Infinite Busbar Model




The single 3-rd order machine infinite busbar model can be established by connecting the 3rd-order machine module in

FIG. 4

with the tieline infinite busbar network module (

FIG. 13

) [25]. The reference frame is stationary, that is Δ, Ω


ref


=0.



Claims
  • 1. A model for testing an entire power transmission system under small perturbation stability comprising:a microprocessor for conducting stability analyses of said power system; a plurality of first type of plug-in transfer blocks representing various components of said power system and processable by said microprocessor wherein each block is derived from sets of equations representing one of said components; a recombination network represented by a second type of transfer blocks, wherein said second type of block is derived from sets of network equations representing the power system network and said second type of block is adapted to receive a first set of said first type of plug-in blocks to form a model of a power system, such that an entire power system may be tested for small perturbation stability by simple matrix manipulation when processed by said microprocessor.
  • 2. A method for modeling a power transmission system under small perturbation stability comprising the steps of:a. providing a microprocessor for conducting stability analyses of said power system; b. preparing a plurality of first type of plug-in transfer blocks representing various components of said power system and processable by said microprocessor wherein each block is derived from sets of equations representing one of said components; c. providing a recombination network represented by a second type of transfer blocks, wherein said second type of block is derived from sets of network equations representing the power system network and said second type of block is adapted to receive a first set of said first type of plug-in blocks to form a model of a power system, such that an entire power system may be tested for small perturbation stability by simple matrix manipulation d. matrix manipulating of said recombination network to test for small perturbation stability of an entire power system.
US Referenced Citations (2)
Number Name Date Kind
5594659 Schlueter Jan 1997
5745368 Ejebe et al. Apr 1998
Non-Patent Literature Citations (10)
Entry
Parniani et al., “Computer Analysis of Small-Signal Stability of Power Systems Including Network Dynamics”, IEE Proceedings Gen., Trans. and Dist., vol. 142, Issue: 6, pp. 613-617, Nov. 1995.
Sanchez-Gasca et al. “Identifying Linear Models from Time Domain Simulations”, IEEE Computer Applications in Power, vol. 10, Issue: 2, pp. 26-30, Apr. 1997.
Yuan et al., “Eigen Analysis Method for Stability Evaluation of Linear Constant-Coeffient RDDE and its Application to Power System”, Proc. IEEE Reg. 10 Conf. on Computer, Comm., Control and Power Eng., vol. 5, pp. 155-156, Oct. 1993.
Wong et al, “Eigenvalue Analysis of Very Large Power Systems”, IEEE Transactions on Power Systems, vol. 3, Issue: 2, pp. 472-480, May 1988.
Campagnolo et al., “An Efficient and Robust Eigenvalue Method for Small-Signal Stability Assessment in Parallel Computers”, IEEE Transactions on Power Systems, vol. 10, Issue:1, pp. 506-511, Feb. 1995.
Angelidis et al., “Efficient Calculation of Critical Eigenvalue Clusters in the Small Signal Stability Analysis of Large Power Systems”, IEEE Transactions on Power Systems, vol. 10, Issue:1, pp. 427-432, Feb. 1995.
Hakavik et al., “Power System Modelling and Sparse matrix Operations Using Object-Oriented Programming”, vol. 9, Issue:2, pp. 1045-1051, May 1994.
Parniani, M.; Iravani, M.R.;“Computer Analysis of Small-signal Stability of Power Including Network Dynamics”, IEE Proceedings—Generasion, Transmission and Distribution, Vo. 142, Issue 6, pp. 613-617, Nov. 1995.
Sanchez-Gasca, J.J.; Clark, K.; Miller, N.W.; Okamoto, H.; Kurita, A.; Chow, J.H.; “Identifying Linear Models from Time Domain Simulations”, IEEE Computer Applications in Power, vol. 10, Issue 2, pp. 26-30, Apr. 1997.
Ben-Tao Yuan; Dao-Zhi Xia; Yi-Xin Ni; “Eigen Analysis Method for Stability Evaluation of Linear Constant-coefficient RDDE and its Application to Power System”, Proceedings TENCON '93, IEEE Region 10 Conference on Computer, Communications, Control and Powe, Oct. 1993.