Stability of Distribution Systems With A Large Penetration of Distributed Generation

Stability of Distribution Systems with a Large Penetration
of Distributed Generation
der Fakultät für Elektrotechnik und Informationstechnik
der
Universität Dortmund
genehmigte
DISSERTATION
zur Erlangung des akademischen Grades
Doktor der Ingenieurwissenschaften
von
Bin Huang M.Sc.(Eng)
Dortmund
Referent: Prof. Dr.-Ing. E. Handschin
Korreferent: Prof. Dr.-Ing. Dr.-Ing. S. Kulig
Tag der mündlichen Prüfung: 9. November 2006

Acknowledgements
My sincere gratitude and appreciation go to my supervisor, Prof. Dr.-Ing. Edmund Hand-
schin, for his skilled guidance, valuable comments, and encouragement that he gave me
over the past years.
Special thanks go to my Korreferent, Prof. Dr.-Ing. Dr.-Ing. Stefan Kulig, whose great
help has been indispensable for fulfilling this work.
I am grateful to all the members of the institute for contributing to such an inspiring
and pleasant atmosphere.
The financial support from the Graduate School of Production Engineering and Logis-
tics at the University of Dortmund is gratefully acknowledged.
Finally, I wish to express my utmost gratitude to my wife, my parents, my sisters, my

brothers-in-law, and my little nephew, whose love, belief and encouragement make me
able to achieve my goal.
i
Contents
List of the Acronyms iii
1 Introduction 1
1.1 Motivation and Objectives of this Work . . . . . . . . . . . . . . . . . . 1
1.2 Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2 Concepts, Attributes and Technologies of Distributed Generation and Power

System Stability 5
2.1 Distributed Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Concept and Attributes of Distributed Generation . . . . . . . . . 5
2.1.2 Overviews of Distributed Generation Technology . . . . . . . . . 6
2.1.3 Impact of Distributed Generation to Power System Operation . . 7
2.2 Power System Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.1 Overview of Power System Dynamics . . . . . . . . . . . . . . . 8
2.2.2 Concepts and Methods for Analyses of Power System Stability . . 9
2.2.3 Simulation Tools for Analyses of Power System Stability . . . . . 13
3 Study of the Eigenvalues of the State Space Model of Distribution Networks 14

3.1 Modeling of Lines/Cables . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.1.1 Selection of a Common Reference Frame . . . . . . . . . . . . . 15
3.1.2 Relation of Network Eigenvalues under abc and dq0 Frames . . . 16
3.1.3 Modeling Lines/Cables in dq0 Reference Frame . . . . . . . . . 20
3.2 Forming of the State Matrix of Distribution Networks . . . . . . . . . . . 23
3.3 A Brief View of the Study Networks . . . . . . . . . . . . . . . . . . . . 26
3.4 Simulations and Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
ii
3.4.1 Properties of Network Eigenvalues . . . . . . . . . . . . . . . . . 27

3.4.2 Analysis of the Characteristics of Distribution Network Dynamics 36
3.4.3 Conclusions and Comments . . . . . . . . . . . . . . . . . . . . 36
4 Impacts of DG on the Small Signal Stability of Distribution Systems 39

4.1 Implementation of the Small Signal Stability Analysis in Power Systems . 40
4.1.1 Basic Concepts of Modal Analysis . . . . . . . . . . . . . . . . . 40
4.1.2 Selection of a Common Reference Frame . . . . . . . . . . . . . 42
4.1.3 Some Specifics of the Power System Model . . . . . . . . . . . . 42
4.2 Models of Power System Components . . . . . . . . . . . . . . . . . . . 43
4.2.1 Equivalent of the Transmission System . . . . . . . . . . . . . . 43
4.2.2 Network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2.3 Constant Speed Wind Turbine Generation System . . . . . . . . . 44
4.2.4 STATCOM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2.5 Loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3.1 Impacts of WGS on the Small Signal Stability of the SNI . . . . . 53
4.3.2 Strategies to Increase the Penetration Level of DG and the Corre-
sponding Impacts on the Small Signal Stability of the SNI . . . . 65
4.3.3 Relevant Results for the SNII . . . . . . . . . . . . . . . . . . . 70
4.3.4 Conclusions and Comments . . . . . . . . . . . . . . . . . . . . 72
5 Impacts of DG on the Voltage Stability of Distribution Systems 74

5.1 Implementation of the Continuation Method . . . . . . . . . . . . . . . . 75
5.2 Models of Power System Components . . . . . . . . . . . . . . . . . . . 80
6 Conclusions and Future Works 86

6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
6.2 Future Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
Appendix 89
Bibliography 100
iii
List of the Acronyms

ASPAS Adaptive Set Point Adjustment Scheme
CCP Common Coupling Point
CE Characteristic Eigenvalues
CN Connecting Node
CRE Critical Eigenvalues
CRF Common Reference Frame
DAE Differential and Algebraic Equations
DG Distributed Generation
EN Ending Node
ODE Ordinary and Differential Equations
OLTC On Load Tap Changer
PWM Pulse Width Modulation
SN Starting Node
SNI Study Network I
SNII Study Network II
TN Trivial Node
UN Upstream Node
WGS Wind Generation System
Chapter 1
Introduction
1.1 Motivation and Objectives of this Work
The large increase of distributed generation (DG), which usually means a small-scale
power generation located at or near to the load being served, has a large impact on system
operations. The growing trend of DG is easily to be recognized because of various politi-
cal, economic and environmental reasons. In order to be prepared for problems that may
happen when the expectation of a high penetration level of DG comes true, huge amount
of efforts of power system engineers have been spent in this field.
It is by no means easy to give a general view of the state of the art for the research
about DG, because of the prompt evolvement and the various emphases of different re-
searches. In the following, a review is presented by focusing on the works which concern
the determination of the allowable penetration level of DG and the impact of DG on the
power system stability.
In [1] the authors attempt to determine a wind power penetration constraint based on
a worst case wind generation change due to a thunderstorm. The authors of [2] dis-
cuss the allowable penetration level of DG from the viewpoint of harmonic distortion
prevention, while the authors of [3] study the similar problem by considering another
factor – the influence of different excitation system control modes on the integration of
distributed synchronous generators. Papers [4][5][6] present the impact of DG on trans-
2 Chapter 1
mission system transient stability, meanwhile, [4][6] also address the transmission system
small signal stability. In papers [7][8][9], the influences of DG on system transient sta-
bility, voltage stability and voltage quality are discussed. Papers [10][11][12][13] focus
mainly on the aspect of islanding operation. In particular, [13] investigate the impact of
high wind power penetration on the frequency stability of power systems under islanding
conditions. In [14], the selection of switch placement to improve the reliability of radial
distribution systems with DG under fault conditions is presented.
As far as the author of this thesis known, to present, the impact of DG on the small
signal stability of distribution systems is relatively less to be discussed. From the author’s
personal viewpoint, the main reasons may be:
• Possibilities of unsymmetric/unbalanced operations of distribution systems -

Compared to transmission systems, distribution systems are subject to higher po-
tential of unsymmetric or unbalanced operations. However, the effective method of
modal analysis for small signal stability mainly reflects the dynamic characteristics
of power systems at a specific equilibrium point under the condition of three phase
symmetry and balanced operations. These constitute a dilemma.
• Lack of a general model of DG - Presently there are several technologies for DG

applications, such as wind power plant, solar photovoltaic generation, fuel cells,
etc. However, even for a certain technology, similar products from different manu-
facturers vary quite a lot about their structure and parameters. This makes the task
of building a general model of DG a demanding challenge. On the contrary, the
small signal stability of the system to a large extent is determined by the eigenval-
ues of the corresponding linearized model, which heavily depend on the structure
and parameter of the system itself. These constitute another dilemma.
Nevertheless, the question has to be answered: to what penetration level of DG can a

distribution system be operated without instability problems?
Thus, one of the objectives of this dissertation is to identify if a high penetration level
of DG will lead to instability problems in a distribution system. If so, the further step is
Chapter 1 3
proposed as to find a solution to eliminate the resulted instability problems; if not, the
further step is proposed as to find some feasible methods to increase the penetration level
of DG.
Power system stability mainly includes categories of transient stability, small signal sta-
bility, voltage stability and frequency stability. Besides small signal stability, this thesis
also gives some discussions about voltage stability. Others, such as transient stability con-
cerning DG and distribution systems, since it has been discussed in many publications,
thus will not be addressed here. As for frequency stability, because it will be not very
likely to happen just due to the integration of many DG in a rather long time period, thus
will not be studied in this work.
When discussing the small signal stability, the following assumptions are adopted:
1. System parameters are three phase symmetric and system operations are three phase
balanced. This may be acceptable at least for German grid, since its distribution
systems are well planned and structured.
2. Among lots of DG technologies, constant speed wind generation systems are specif-
ically designated to proceed the investigation, since, on the one side, they have
a relatively strong interaction with the connected system; on the other side, their
mathematically mature model is available.
Another objective of this thesis is to give an answer to the following consideration: Com-
pared to transmission system, disturbances occur in distribution system more often. The
connection of DG will make this situation even worse. Hence a reasonable question can
be as follows: will these apparently continuous disturbances act as a chain to cause some
dangerous events in the electrical network itself (dynamics of other components, such as
generators and DG, are not included)?
This question is answered by investigating the relation between the eigenvalues of the
state space model of the electrical network and the reasonable expectations of the time
interval of two successive disturbances.
4 Chapter 1
The obtained results can simultaneously serve as a verification and validation of the cor-
rectness of the network model: As well known, when discussing power system stability,
normally the network is mathematically modeled as a constant nodal admittance matrix,
i.e., the network is treated in steady state, because its dynamics die out quite fast. Now,
facing to a distribution system, since the parameters of the lines/cables differ a lot from
those of a transmission system, the assumption is still correct, i.e., the network dynamics
are still negligible?
1.2 Outline of the Thesis

Following the introduction given in chapter 1, chapter 2 overviews the technical back-
ground of this work, mainly including:
1. the concept and attributes of DG, the state of the art of their technologies and their
impacts on system operations.
2. various aspects of power system stability, such as the classification of them, the
applicable analysis methods and the available simulation tools.
Chapter 3 focuses solely on the dynamics of the electrical network (without other dynamic
components, such as synchronous generators, induction motors, etc). The whole network
model is built under dq0 reference frame after a discussion of the relations between the
abc and dq0 coordinations. Consequently the network dynamics are studied by evaluating
the eigenvalues of the corresponding state matrix.
Chapter 4 investigates the small signal stability of the whole system based on modal
analysis when the penetration level of DG achieves a high extent.
Chapter 5 concentrates on the influence of DG and its operation modes on the voltage
stability. During the study, the continuation power flow method is adopted to determine
the system loadability under different scenarios.
Finally, chapter 6 summarizes the whole work and points out the potential directions for
future research.
Chapter 2
Concepts, Attributes and Technologies

of Distributed Generation and Power
System Stability
2.1 Distributed Generation
2.1.1 Concept and Attributes of Distributed Generation
During the past several decades, power generation has been highly centralized in large
power plants. Customers are served primarily by utility distribution companies that have
connections to large generation facilities using high-voltage transmission lines and con-
nections to customers through their lower-voltage distribution lines.
This operating structure is built on the basis of economy, security, and quality of sup-
ply. It is operated by hierarchical control centers and allows the system to be monitored
and controlled continuously. The generation is instantly adjusted to the consumption by
monitoring the frequency and on the basis of elaborated load forecasting models. The
voltage is controlled to be within specific limits by means of appropriate coordinated de-
vices, generators, on load tap changers (OLTC), reactive compensation devices, etc.
Nowadays there is a considerable interest in connecting generation to the distribution

system and this is termed as DG.
6 Chapter 2
Compared to conventional generation, DG usually has characteristics as [15][16]:
• not centrally planned

• not centrally dispatched
• with small capacity
• usually connected to the distribution system.
And it mainly contributes to [15][16]:
• reduction of emission levels (mainly CO2)

• energy efficiency or rational use of energy
• deregulation or competition policy
• diversification of energy sources
• ease of finding sites for smaller generators
• short construction times and lower capital costs of smaller plant
• reduction of power transmission costs and losses.
2.1.2 Overviews of Distributed Generation Technology

The currently available DG technologies up to a capacity of 5 MW, mainly include [17]:
• Microturbines - This new and emerging technology is currently only available

from a few manufacturers. Other manufacturers are looking to enter this emerging
market, with models ranging from 30 to 200 kW. Microturbines promise low emis-
sion levels, but the units are currently relatively expensive. Obtaining reasonable
costs and demonstrating reliability will be major hurdles for manufacturers. Micro-
turbines are just entering the marketplace, and most installations are for the purpose
of testing the technology.
• Industrial combustion turbines - This mature technology ranges from 1 MW to

over 5 MW. They have low capital cost, low emission levels, but also usually low
electric efficiency ratings. Development efforts are focused on increasing efficiency
levels for this widely available technology. Industrial combustion turbines are being
used primarily for peaking power and in cogeneration applications.
Chapter 2 7
• Fuel cells - Although the first fuel cell was developed more than 150 years ago, this
technology remains in the development stage. Fuel cell emission levels are quite
low, but cost and demonstrated reliability remain major problems for its market
penetration. The few fuel cells currently being used provide premium power or are
in applications subsidized by the government or gas utilities.
• Photovoltaics - Commonly known as solar panels, photovoltaic (PV) panels are

widely available for both commercial and domestic use. Panels range from less
than 5 kW and up to several MW. They produce no emissions, and require minimal
maintenance. However, they can be quite costly. Less expensive components and
advancements in the manufacturing process are required to eliminate the economic
barriers now impeding the wide-spread use of PV systems. Photovoltaics may be
used in remote locations without grid connections.
• Wind turbine systems - Wind turbines are currently available from many man-
ufacturers. Compared to other renewable technologies, they provide a relatively
inexpensive way to produce electricity. But as relying upon the variable and some-
what unpredictable wind, they are not suited for continuous power production.
2.1.3 Impact of Distributed Generation to Power System Operation

From a technical aspect, DG poses some impacts on power systems, especially to the
operations of distribution systems. They mainly include [15][16][18]:
• Voltage Issue - The utility must provide electricity to the customers at a voltage
within specific limits. However, if the capacity of the connected DG is relatively
large or the connection between transmission and distribution system is relatively
weak, the steady state voltage rise may be a problem.
• Protection - The connection of DG changes the magnitude, duration, and direction

of a fault current, thus requires corresponding adaptations of the existing protection
system.
• Power Quality - Most commonly recognized power quality events include fre-
quency fluctuation, sustained over voltages, high frequency noise, harmonic distor-
tions, high voltage spikes and voltage sags and swells. DG may decrease the power
8 Chapter 2
quality when it has relatively large capacity or is connected to a weak distribution

system.
• Stability - The increment of the capacity of DG determines that it plays a more and
more important role in power system stability. Therefore more considerations are
required in this aspect.
• Security - DG energizes the system via lots of points thus complicates the policies
of isolation and earthing for safety before maintenance work is undertaken.
2.2 Power System Stability
2.2.1 Overview of Power System Dynamics

In order to gain a better understanding of power system stability, some knowledge of the
dynamic phenomena in power system is necessary. Dynamics here means the character-
istics of the way the power system responds to a disturbance. Such a disturbance can be
either an intentional operation, e.g., a scheduled generator switching, or an accident, e.g.,
a light striking. One of the common and perhaps most important classification of power
system dynamics is their natural time range of response. There are various categorizations
of power system dynamic phenomena by considering them in different levels of detail. In
this work, the classification given in [19] is adopted.
Power system dynamics can be divided into four groups as shown in fig. 2.1. Among
them, wave dynamics, or surges, are the fastest dynamics, which occur in high volt-
age transmission lines and result from the propagation of electromagnetic waves caused
by lightning strikes or switching operations. Slower are the electromagnetic dynamics,
which take place in the machine windings following a disturbance. Much slower are the
electromechanical dynamics, which describe the oscillation of the rotating masses of gen-
erators and motors when the system is subject to disturbances. Thermodynamics are the
slowest ones which correspond to phenomena such as boiler control action in steam power
plants as a response to frequency variations. In the following sections, the discussed
problems of power system stability mainly belong to the categories of electromechanical
dynamics and thermodynamics, i.e., they vary from seconds to many minutes.
Chapter 2 9
Figure 2.1: Time range of power system dynamics [19]
2.2.2 Concepts and Methods for Analyses of Power System Stability
Stability, as a fundamental concept in most engineering disciplines, has been extensively

studied and more than 28 definitions of it are introduced in the system theory focused
on different technical and physical aspects. In power system engineering, the proposed
definition, in a recent work of IEEE/CIGRE Joint Task Force [20], is formed as:
• Power system stability is the ability of an electric power system, for a given initial
operating condition, to regain a state of operating equilibrium after being subjected
to a physical disturbance, with most system variables bounded so that practically
the entire system remains intact.
In the early stage, power system engineers mainly encountered instability problems re-
lated to generator rotor angle stability, with which the whole system can not keep syn-
chronism after a disturbance. However, with the continuous expansion of modern power
systems under more stressed operation conditions, different forms of system instability
have emerged.
A more detailed classification according to the different mechanism is helpful to under-

stand the problems well, though the power system stability is essentially a single problem.
In [20], three categories have been proposed, i.e., rotor angle stability, voltage stability and
frequency stability.
10 Chapter 2
Rotor Angle Stability
Rotor angle stability refers to the ability of synchronous machines of an interconnected

power system to remain in synchronism after being subjected to a disturbance. It depends
on the ability to maintain/restore equilibrium between electromagnetic torque and me-
chanical torque of each synchronous machine in the system. Instability that may result
occurs in the form of increasing angular swings of some generators leading to their loss
of synchronism with other generators [20].
By considering the size of the disturbance, it can be further divided into two subcate-
gories:
• Small disturbance rotor angle stability, sometimes also referred to as small signal
stability, is concerned with the ability of the power system to maintain synchronism
under small disturbances, such as a small load increment. Small signal stability
depends on the initial state of the system. In today’s power systems, small signal
instability problems are usually caused by the lack of enough damping torque, un-
der which some variables of the power system show an oscillation characteristics
with increasing amplitude after being subjected to a disturbance.
Since the disturbance is assumed to be small, the linearization of system model

is allowable and thus the modal analysis methods are applicable [21][22][23]. Non-
linear time domain simulation is another effective approach in the analysis of the
small signal stability problems. The results of modal analysis usually have to be
verified by using the time domain simulations as a benchmark. However, the sim-
ulation method has relatively weak ability to determine the critical modal charac-
teristics, which is important for designing appropriate countermeasures; while the
modal analysis gives complete information on the inherent characteristics of modes
(mode shapes, participation factors, transfer functions, etc.), which can be used in
control design for damping enhancement. A feasible proposal is to use both of them
in a complementary manner [24]. Besides, bifurcation theory is applicable by notic-
ing that the small signal stability is bounded by the occurrence of Hopf bifurcation
(where a pair of complex-conjugate eigenvalues of the system appears) [25][26].
Chapter 2 11
• Large disturbance rotor angle stability, sometimes also named as transient stability,
is concerned with the ability of the power system to maintain synchronism when
subjected to a severe disturbance, such as a short circuit on a transmission line.
Transient stability depends on both the initial operating state of the system as well
as the severity of the disturbance. Instability problems are usually due to insuffi-
cient synchronizing torque thus evolve normally in the form of aperiodic angular
deviations.
The highly nonlinear characteristics of transient stability decide that the nonlin-
ear time domain simulation will play an important role in relevant studies. This
simulation is actually a numerical procedure to solve the differential and algebraic
equations (DAE) of a power system [21][27]. Another way to determine the stabil-
ity is the so called direct method, in which a Lyapunov function is formed and thus
Lyapunov’s second method is adopted for the evaluation. However, this method
mainly has two difficulties: 1. How to construct a good Lyapunov function; 2. How
to assess a practical stability domain [21][27][28]. Further, methods based on pat-
tern recognition are also proposed for transient stability analysis, which try to infer
results about the present stability properties by digging information from the ”past
experience” [27].
Voltage Stability
Voltage stability refers to the ability of a power system to maintain steady voltages at all
buses in the system after being subjected to a disturbance from a given initial operating
condition. It depends on the ability to maintain/restore equilibrium between load demand
and load supply from the power system. Instability that may results occurs in the form
of a progressive fall or rise of voltages of some buses. A possible outcome of voltage
instability is loss of load in an area, or tripping of transmission lines and other elements
by their protective systems leading to cascading outages [20].
The loads are usually the driving force for voltage instability. When subjected to a dis-
turbance, the power consumption of the loads tends to be restored due to the action of
OLTC, thermostats, etc. This leads to the reactive power consumption increment thus
12 Chapter 2
further voltage reduction. A progressive drop of the voltage persists which finally results
in voltage collapse, when there is no enough generator capability or network transfer ca-
pacity available. The network transfer capacity is determined by the voltage level and the
magnitude of the corresponding impedance [21][29].
As in transient stability analysis, time domain simulation also serves as an effective and
accurate tool for analyzing voltage stability. However the computation requirement is
demanding. In order to speed the analysis, some simplifications are made by assuming
the system in a steady state thus those static methods are applicable, such as using con-
tinuation method to determine the loadability limits and adopting some sensitivity index
to estimate the distance from the present system operation point to the voltage collapse
point [21][29][30]. Bifurcation analysis is also appropriate for voltage stability since it
is observed that the instability that may occur usually coincides with the singularity of
system state Jacobian, which is the necessary condition of saddle-node bifurcation [29].
Frequency Stability
Frequency stability refers to the ability of a power system to maintain steady frequency
following a severe system upset resulting in a significant imbalance between generation
and load. It depends on the ability to maintain/restore equilibrium between system gen-
eration and load, with minimum unintentional loss of load. Instability that may result
occurs in the form of sustained frequency swings leading to tripping of generating units
and/or loads [20].
The possible instability is usually analyzed with time domain simulation. Though the
procedure is similar as that in transient stability simulation, the model of the power sys-
tem may need some adaptation. For example, since the potential frequency instability is
usually a long term dynamic process after the system being subjected to a severe upset,
proper representations of the prime mover, energy supply systems and functions of wide
range protection and control systems may be necessary [21].
Chapter 2 13
Comments on Classification
The classification of power system stability helps identifying causes of instability, ap-
plying suitable analysis tools, and developing corrective measures. However, in a real
situation any one form of instability may not occur solely. This is true especially for
highly stressed systems and for cascading events. Thus power system stability will by all
means be considered as a whole; no category of the mentioned stability would be bettered
at the expense of another [20].
2.2.3 Simulation Tools for Analyses of Power System Stability

Nowadays many powerful software packages are available for analyses of power system
stability. Some of them are mature and commercial tools, while the others are useful and
open source programs.
PSAT, a freely distributed Matlab-based power system analysis program developed by Dr.
Federico Milano [31], stands out among those open source and freely available programs
for its relatively complete functionalities. It includes functions as power flow, continu-
ation power flow, optimal power flow, small-signal stability analysis, and time-domain
simulation. To fulfill these analyses, PSAT supports a variety of static and dynamic mod-
els. For example, in a static analysis, a synchronous generator can be modeled as a slack
bus or a PV bus; however, in a dynamic analysis, options of the synchronous generator
model with a dynamic order varying from 2 to 8 are provided. Besides, PAST includes a
user friendly graphical interface and a Simulink-based one-line network editor [32].
These characteristics make PSAT attractive thus be adopted for the work presented in this
thesis. During the research, all concerned codes have been carefully investigated. Some
small errors have been fixed and some models have been adapted/developed. About these
aspects, more details are given in following chapters.
Chapter 3
Study of the Eigenvalues of the State

Space Model of Distribution Networks
Network dynamics are normally omitted in order to simplify power system stability anal-
ysis. This is based on the observation that the components of power systems have quite
different time constant of the natural response, for which a general view is given in [33].
[34] further justifies the elimination of network transients in power system stability anal-
ysis via using singular perturbation theory. Though papers [35][36][37] put forward sug-
gestions and accompanying reasons of including network dynamics into a large scale
analysis, considering this will prominently increase the dynamic order of the studied sys-
tem, the feasibility of it is still open for more investigations.
The necessity of including network dynamics is studied in this thesis based on the follow-
ing observations: firstly, the trend of integrating more and more DG into the present power
systems is well recognized; secondly, compared to traditional electricity sources/loads
connected to the distribution network, these DG have a relatively large capacity up to
several megawatt; last but not least, the operation of DG has more or less some stochastic
properties, for which the output of wind generation systems or photovoltaic systems may
serve as an example. All these combined together, leads to the question studied here, i.e.,
if the incessant stochastic inputs (power, currents ... etc.) will act as a chain to induce
harmful responses (nodal overvoltage ... etc.) in the distribution network? This question
has been investigated in this chapter by analyzing the properties of the eigenvalues of the
distribution network state space model.
Chapter 3 15
In this chapter, section 3.1 focuses on modeling lines/cables in dq0 reference frame. Sec-
tion 3.2 introduces a method to form the state space model of a distribution network.
Information about the study networks is given in 3.3, while section 3.4 presents simula-
tions, results and comments.
The reason of modeling the network in dq0 reference frame is, if the analysis results
demonstrate that the negligence of the distribution network dynamics is inappropriate,
the developed models and programs can be directly incorporated to other system compo-
nent models for the consequent analysis.
Main conclusions of this chapter are:
1. For distribution networks, since generally the dynamics die out quite fast, the ap-
parently successive stochastic disturbances have nearly independent influences on
the network. Hence these impacts can be studied independently.
2. It is not necessary to include the electric network dynamics into the stability anal-
ysis. The distribution network can still be mathematically modeled via a nodal
admittance matrix.
3.1 Modeling of Lines/Cables
3.1.1 Selection of a Common Reference Frame

The introduction of Park Transformation simplifies the mathematical description of syn-
chronous machines. It defines a new set of stator variables such as currents and voltages,
which are obtained by projecting the actual variables on abc axes to dq0 axes. However,
during the analysis, each synchronous machine is modeled in its own dq0 reference frame.
When simulating an interconnected system, all the variables have to be referred to a com-
mon reference frame (CRF). After the determination of the CRF, other system variables
and component models, e.g., models of induction motors, normally are also referred to
this frame to simplify the analysis [21][38][39].
16 Chapter 3
To build the model of the electrical network under dq0 frame, Kron transformation [40]
is adopted. In this work, the reference frame of Kron transformation is selected the same
as the CRF which is specifically designated as: in which, the q axis leads the d axis by
90◦ and coincides with the system reference phasor, as shown in fig. 3.1.
Figure 3.1: Selection of the common reference frame
Hence the transformation matrix from abc reference frame to dq0 reference frame is
 
sin θ sin(θ − 2π3 ) sin(θ + 2π
3 )
2 2π 2π

P =  cos θ cos(θ − 3 ) cos(θ + 3 ) 
 (3.1)
3 
1 1 1
2 2 2
where
θ = ω0t, ω0 = 2π f , f the system frequency.
The inverse transformation matrix is

 
sin θ cos θ 1
P−1 = 
 
2π 2π (3.2)
 sin(θ − 3 ) cos(θ − 3 ) 1 

sin(θ + 2π 2π
3 ) cos(θ + 3 ) 1
3.1.2 Relation of Network Eigenvalues under abc and dq0 Frames

To investigate the properties of network dynamics by evaluating the eigenvalues of the
network under dq0 frame, the theoretical precondition is that the eigenvalues under abc
frame can be reconstructed from those under dq0 frame. This can be achieved based on
the following observation:
Chapter 3 17
Observation
Assuming a symmetric network has 3N state variables, thus it has 3N eigenvalues. Under
dq0 reference frame, these 3N eigenvalues consist of N pairs of eigenvalues which are
complex-conjugate/double-real ones and relate only to the dq components (lately referred
to as λdq ), and N eigenvalues which relate only to the 0 components (lately referred to as
λ0 ). Correspondingly, under abc frame the 3N eigenvalues can also be divided into two
groups, one of which contains N pairs of eigenvalues (lately referred to as λabc2N ), the
other includes the rest (lately referred to as λabcN ). Then the eigenvalues of the network
model in abc and dq0 reference frame have the following relations:
• For λabc2N , assuming it includes complex-conjugate eigenvalues σ1 ± jω1 , after

dq0 transformation, the corresponding complex-conjugate eigenvalues will be σ1 ±
j(ω1 − ω0 ) and σ1 ± j(ω1 + ω0 ).
• For λabcN , this set of eigenvalues will keep unchanged after dq0 transformation.
Explanation
Reference [40] describes the relations of the eigenvalues of electrical networks under abc
and dq0 frames. But the discussion is just confined to the basic network components such
as a series impedance and a shunt capacitance. A general conclusion is not available if the
relations presented above are applicable to a much more general and complicate network,
for example, the one given in fig. 3.6.
It can be foreseen that a strictly mathematic proof may be too complicated and ardu-
ous to achieve. In this thesis, this question is discussed from another aspect, which is
much more based on a physical meaning rather than a mathematical reasoning.
The discussion consists of two steps.
• Step 1
The set of λabc2N is assumed to have complex-conjugate eigenvalues σ1 ± jω1 (from
the later example it can be seen these complex-conjugate eigenvalues usually ap-
pear as two pairs).
18 Chapter 3
Firstly, the input signals Iabc in fig. 3.2 are assumed to be positive sequential si-
nusoids with frequency of ω1 . It is well known that when the input signals are
applied to the network at t = 0, the output signals Oabc can be considered as a com-
bination of transient signals as well as steady signals. The transient signals will
die out as time is passing, while the steady signals will exist as long as the input
signals sustain. The amplitude of the steady signals are determined by both the
poles and the zeros of the system. For simplifying the discussion, it is assumed the
steady output signals have a prominently large amplitude, i.e., there exists a res-
onance in the responses. It can be shown that after dq0 transformation, the input
signals Idq0 in fig. 3.2 consist of only dq components, which are sinusoidal signals
with frequency of ω1 − ω0 , while the 0 component constantly equals to zero [41].
Meanwhile, when applying dq0 transformation to the outputs signals Oabc , which
is equivalent to apply the transformation separately to the transient signals as well
as the steady signals in Oabc , this leads to: for the transient signals, their decay-
ing ratio will keep unchanged, which means the real parts of the eigenvalues hold
unchanged; for the steady signals, the frequency will be shifted to ω1 − ω0 , which
means under dq0 frame the resonance frequency is shifted from ω1 to ω1 − ω0 .
Combining these together indicates that the network model under dq0 frame has
eigenvalues at σ1 ± j(ω1 − ω0 ).
Then negative sequential sinusoids with frequency of ω1 are assumed to be the

input signals Iabc . It is noted that the negative sequence signals have nearly the
same impact to the network as that of the positive sequence signals. Following a
similar reasoning as mentioned above, it can be found that the network model under
dq0 frame has eigenvalues at σ1 ± j(ω1 + ω0 ).
It is worth noting that the dq0 transformation will also alter the frequency of the
transient signals. However, in above discussions the attention is only focused on
the alteration of the decaying ratio of the transient signals. This may help to make
the discussions easier.
• Step 2
Chapter 3 19
Suppose that the set of λabcN has complex conjugate eigenvalues σ2 ± jω2 and
the input signals Iabc are three phase identical (zero sequence) and sinusoidal with
frequency of ω2 . A similar deduction procedure as in step 1 reveals that the eigen-
values hold the same under abc and dq0 frames.
It may help to note that during the deduction, the application of dq0 transforma-
tion to three phase identically sinusoidal signals will lead to zero in dq axis while a
sinusoidal signal the same as the input signal in 0 axis.
Figure 3.2: System under abc and dq0 frames
Example
A small network shown in fig. 3.3 is used to illustrate the given observation. Its parame-
ters are as follows:
Impedance
Node Node Rs Rm Ls Lm
1 2 0.02 0.012 0.2 0.1
2 3 0.04 0.01 0.35 0.18
3 4 0.08 0.035 0.27 0.13
3 5 0.1 0.06 0.5 0.23
Capacitance
Node Cs Cm
2 0.01 0.002
The eigenvalues under abc and dq0 frames are given in table 3.1.
20 Chapter 3
Figure 3.3: Example for investigating the network eigenvalues under abc and dq0 frames
λabc2N −18.058 ± j13055 −18.058 ± j13055 −73.708

λabc −73.708 −45.355 −45.355
λabcN −18.697 ± j4932.4 −78.436 −41.588
λdq −18.058 ± j13369 −18.058 ± j12740 −73.708 ± j314.16
λdq0 −45.355 ± j314.16
λ0 −18.697 ± j4932.4 −78.436 −41.588
Table 3.1: Eigenvalues of the system given in fig. 3.3
3.1.3 Modeling Lines/Cables in dq0 Reference Frame
Before deducing the network dynamic model, the lines/cables are modeled in advance.
A single phase equivalent of a serial impedance is shown in fig. 3.4(a). In abc frame,
it is described by
dIabc
L + RIabc = V1abc − V2abc (3.3)
dt
where Iabc = [ia ib ic ]T , Vabc = [va vb vc ]T . The superscript T means transpose. Besides,
since the network is assumed to be symmetric, thus L and R are
   
Ls Lm Lm Rs Rm Rm
   
L=
 Lm Ls Lm  , R =  Rm Rs Rm 
  
Lm Lm Ls Rm Rm Rs
where
Ls is the self inductance, Lm the mutual inductance
Rs the self resistance, Rm the mutual resistance
Similarly, a single phase equivalent of a shunt capacitance is shown in fig. 3.4(b). It

Chapter 3 21
is described by
dVabc
C = I1abc − I2abc (3.4)
dt
with  
Cs Cm Cm
 
C=
 Cm Cs Cm


Cm Cm Cs
where Cs is the self capacitance, Cm the mutual capacitance.
Figure 3.4: Single line diagram of basic components of lines/cables
By applying Kron transformation to eqn. (3.3) leads to

d(P−1 Idq0 )
PL + PRP−1 Idq0 = V1dq0 − V2dq0 (3.5)
dt
As discussed in [41], eqn. (3.5) can be rearranged as

did
L − ω0 Liq + Rid = v1d − v2d
dt
diq
L + ω0 Lid + Riq = v1q − v2q
dt
di0
L0 + R0 i0 = v10 − v20 (3.6)
dt
where
L = Ls − Lm , L0 = Ls + 2Lm
R = Rs − Rm , R0 = Rs + 2Rm
By applying Kron transformation to eqn. (3.4) leads to

dP−1 Vdq0
PC = I1dq0 − I2dq0 (3.7)
dt
22 Chapter 3
which can be rearranged as

dvd
C − ω0Cvq = i1d − i2d
dt
diq
C + ω0Cvd = i1q − i2q
dt
dv0
C0 = i10 − i20 (3.8)
dt
where
C = Cs −Cm , C0 = Cs + 2Cm
When choosing the following base quantities for the variables in dq0 frame
vb = peak value of rated line-to-neutral voltage, V
ib = peak value of rated line current, A
fb = rated frequency, Hz
Then the per unit equations of eqn. (3.6) is:
L̄ d īd
− L̄īq + R̄īd = v̄1d − v̄2d
ω0 dt
L̄ d īq
+ L̄īd + R̄īq = v̄1q − v̄2q
ω0 dt
L̄0 d ī0
+ R̄0 ī0 = v̄10 − v̄20 (3.9)
ω0 dt
where the super bar denotes a per unit variable.
Similarly the per unit equation of (3.8) is
C̄ d v̄d
− C̄v̄q = ī1d − ī2d
ω0 dt
C̄ d v̄q
+ C̄v̄d = ī1q − ī2q
ω0 dt
C̄0 d v̄0
= ī10 − ī20 (3.10)
ω0 dt
Since the system is assumed to operate under balanced conditions, the zero sequence volt-
ages and currents in eqn. (3.9) and (3.10) keep to trivial thus the corresponding equations
Chapter 3 23
are dropped. This leads to

L did
− Liq + Rid = v1d − v2d
ω0 dt
L diq
+ Lid + Riq = v1q − v2q (3.11)
ω0 dt
and
C dvd
−Cvq = i1d − i2d
ω0 dt
C diq
+Cvd = i1q − i2q (3.12)
ω0 dt
where the super bar is dropped to make the expressions concise.
3.2 Forming of the State Matrix of Distribution Networks

[42] and [43] have introduced some general algorithms of obtaining the state matrix of a
general circuit. Although they are capable of treating the electrical network of a power
system, the development of such a program is not an easy work. In this thesis, the specif-
ically studied object is distribution networks, which normally are weakly meshed. Based
on this attribute, a special way of forming the state matrix is proposed.
The implementation of the proposed method is explained with the simple network shown
in fig. 3.5. An assumption is adopted here that in the example between two nodes there ex-
ists only an impedance. The capacitance of the lines/cables will be presented and treated
as an independent shunt component.
The procedure consists of six steps as:
Step 1, Nodal Characteristics Analysis

For a distribution network, three kinds of nodes are important for the analysis. They are:
1. Starting Node (SN). SN is usually the common coupling point (CCP) of a transmis-
sion network and a distribution network. In fig. 3.5, node 1 is the SN.
2. Ending Node (EN). EN is a node where only one transmission line is connected. In
fig. 3.5, node 7, 8, 9, 10 and 11 are EN.
24 Chapter 3
Figure 3.5: Example for explanation of the procedure to form network state matrix
3. Connecting Node (CN). CN is a node which connects at least two lines/cables as

well as other energy sources/sinks (such as DG, shunt capacitors, loads, etc). In fig.
3.5, node 2 and 3 are CN. It is worth noting that though node 4, 5 and 6 connect
also more than one transmission lines, they are not CN since there is no other energy
sources/sinks connected to them.
Besides these, for convenience of the subsequent illustration, all remaining nodes will be
referred to as Trivial Node (TN). In fig. 3.5, node 4, 5, 6 are TN.
Step 2, The Selection of State Variables

After the determination of SN, EN CN and TN, the state variables of the network are
specifically designated as:
1. Currents injected into CN and EN. For example, the current flows from node 1 to
node 2 will be designated as a state variable and will be briefly named as i2 rather
than i12 . Then for state variable i3 it denotes the current flowing from node 2 to 3.
Besides, the definition of the positive direction of the current, for energy sources is
that of flowing out, for energy sinks that of flowing in, for transmission lines that of
conforming to a general direction of SN → CN → EN.
Chapter 3 25
2. Voltages of the nodes which have a shunt capacitance branch. In fig. 3.5 only the
voltage of node 2 is chosen as one of the state variables.
Step 3, Nodes Subgrouping

Then the whole network can be divided into subgroups which include the following nodes:
1. CN and EN which are connected upstream via some common TN to a specific CN.
2. The specific CN itself. Afterwards this node will sometimes be referred to as Up-
stream Node (UN).
In fig. 3.5, there exist four subgroups, i.e., (1,2), (2,3), (2,7,8), and (3,9,10,11).
Step 4, Forming a General Form for the Differential Equations of the State Vari-
ables of Branch Currents
For the above mentioned subgroups, as shown in appendix A, each can be described by a
generalized form of differential equations as:
İ2(n−1)×1 = A2(n−1)×2(n−1) I2(n−1)×1 +B2(n−1)×2 v2×1 +C2(n−1)×2(n−1) V2(n−1)×1 (3.13)
where
n is the number of the nodes of respective groups;
I consists of d, q components of the state variables of the currents of the corresponding
group;
v are the d, q components of the nodal voltage of related UN. V contains the d, q compo-
nents of the nodal voltage of that group except the UN.
For example, for group (2,7,8), I is [i7d i7q i8d i8q ]T , v = [v2d v2q ]T , and V = [v7d v7q
v8d v8q ]T .
Step 5, Forming a General Form for the Differential Equations of the State Vari-
ables of Nodal Voltages
For the nodes with a shunt capacitance, they also have a generalized form of differential
equations as:
V̇2×1 = A2×2 V2×1 + B2×2 ∑ I (3.14)
26 Chapter 3
where
V are the d, q components of the nodal voltage;
∑ I is the summation of the currents flowing in/out the node except the one which flows
through the shunt capacitance branch itself.
Taking node 2 as an example, where V is [v2d v2q ]T , and ∑ I equals

 
i2d + i2DGd − i7d − i8d − i3d − i2Loadd
 
i2q + i2DGq − i7q − i8q − i3q − i2Loadq
Step 6, If Necessary, Interfacing the Electrical Network Model with Other System
Components
Equations (3.13) and (3.14) give enough information about the dynamics of the network.
If necessary, a complete system dynamic model can be formed after interfacing them with
the dynamic models of other components, such as generators, DG, etc.
3.3 A Brief View of the Study Networks

In this work, two networks are available for the investigation as:
• Study Network I (SNI) —- which is presented in fig. 3.6. Its structure is based
on IEEE 34 node test feeder [44][45]. Meanwhile, in order to make the network
proper for the study of this thesis, some modifications are adopted as follows
1. The two voltage regulators sited in line segment 814-850 and 852-832 respec-
tively are omitted. The transformer sited in line segment 832-888 is replaced
by a line with the length of 1 km.
2. All lines are assumed to be symmetric.
3. All loads are assumed to be three phase balanced. Distributed loads are sub-
stituted by equivalent spot loads. In general, the magnitudes of the loads are
raised up to about ten times bigger than those in the original case.
More detailed information is given in appendix B.

Chapter 3 27
Figure 3.6: Topology of the SNI
• Study Network II (SNII) —- which is presented in fig. 3.7. The structure and
relevant data are from a real distribution system in Germany. It is worth noting that
this network is connected to the transmission system through both node 1 and node
2. More detailed information is given in appendix C.
Figure 3.7: Topology of the SNII
In the following parts of this thesis, all the proposed methods and analyzing procedures
are explained in detail based on the study of the SNI. However, the obtained conclusions
in the stage of studying the SNI will be verified by the study of the SNII.
3.4 Simulations and Results
3.4.1 Properties of Network Eigenvalues

Before turning to further analysis, the term ’Characteristic Eigenvalues’ (CE) is intro-
duced in advance. Then the distribution of the CE of the SNI and its variants is investi-
gated. Finally the obtained conclusions are verified by studying the distribution of the CE
of the SNII.
28 Chapter 3
Concept of the Characteristic Eigenvalues
Fig. 3.8(a) presents a very simple system: a load is fed by a generator via a transmis-
sion line. For the transmission line, if it has the parameters given in appendix B while
with the length as 1km, then its π equivalent model will have nine eigenvalues under abc
frame as given in fig.3.8(b), which can be separated into two sets, i.e., λabc2N and λabcN .
Since assuming a balanced operation, only the dq circuit components are included to form
the state matrix of the network, thus the corresponding eigenvalues only include λdq , the
counterpart of the set λabc2N . The set of λabc2N can be reconstructed from λdq while λabcN
is dropped, as shown in fig.3.8(c). The CE of the given system are deliberately selected
as in fig.3.8(d). It can be seen that they have included all the necessary information for
the set λabc2N .
Figure 3.8: Explanation for the selection of characteristic eigenvalues
Compared to fig. 3.8(c), it can be noted that in fig. 3.8(d) the pair of zero eigenvalues
is dropped. The reason is they have trivial meaning and their occurrence as well as the
quantities can be expected by investigating the network topology. A further discussion is
given as follows.
Observation
For a network, if its SN, EN and CN are all connected with an equivalent shunt capacitor,
the set λabc2N will consist of two eigenvalues at the origin (the definitions of SN, EN and
CN is in section 3.2).
Discussions
The explanation of the above observation is divided into two steps. Step 1 justifies that if
the requirements are met, the system will have eigenvalues at the origin. Step 2 verifies
Chapter 3 29
the quantity of the zero eigenvalues is 2.
• Step 1:
The state space model of the network is given by ẋ = Ax + Bu. From the analysis
in section 3.2, it can be seen that SN, EN and CN are the interface of the network
to other system components, such as generators, loads, etc. If every SN, EN and
CN has an equivalent shunt capacitance, then in the state space model, the input
vector u only consists of currents (u = [i1 i2 ... in ]T ). To simplify the discussion,
all the currents are set to zero except i1 , i.e., u = [i1 0 ... 0]T (here i1 is a con-
stant value). Under this assumption, since all the system nodes are either isolated
from the ground or connected to the ground via a shunt capacitance, thus theoreti-
cally the constant current source i1 will drive all the nodal voltages increasingly to
infinity. This matches the observation that the network has eigenvalues at the origin.
However, if some nodes of the SN, EN and CN do not have an equivalent shunt
capacitance, then the input vector u will consist of both currents and voltages(u =
[i1 ... vi ... in ]T ). If still specifically chosen u = [i1 0 ... 0]T , this physically means
node i has been directly grounded. Thus the constant current source i1 will only
result in a limited system response (for example, limited nodal voltages). Therefore
there are no eigenvalues at the origin in this case.
• Step 2:
To directly justify the quantity of the zero eigenvalues is difficult. Thus some ex-
amples are studied and the numerical results are reported instead.
Fig. 3.9(a) gives a simple system, in which the network parameters are the same as
those of appendix B. In fig. 3.9(b)-(e), symbols for generators and loads are omitted
to make the plot clear. In fig. 3.9(a), all nodes do not have a shunt capacitance. On
the contrary, all the SN, EN and CN in fig. 3.9(b) have a shunt capacitance. While
in fig. 3.9(c)-(e), one of the SN, EN and CN does not have a shunt capacitance. The
corresponding eigenvalue sets of λabc2N are reported in table 3.2. It can be seen that
only in the case given in fig. 3.9(b), two eigenvalues appear at the origin.
30 Chapter 3
Figure 3.9: Examples for investigating the properties of the characteristic eigenvalues
Case of λabc2N
fig.3.9(a) -80.331, -80.331, -49.716, -49.716, -25.133, -25.133
fig.3.9(b) −37.667 ± j8115.9, −37.667 ± j8115.9, −24.764 ± j6205, −24.764 ± j6205,
−15.159 ± j16598, −15.159 ± j16598, 0, 0
fig.3.9(c) −34.048 ± j7866.8, −34.048 ± j7866.8, −23.51 ± j3972.7, −23.51 ± j3972.7,
−20.032 ± j13372, −20.032 ± j13372
fig.3.9(d) −36.927 ± j7995.1, −36.927 ± j7995.1, −28.096 ± j4706.2, −28.096 ± j4706.2,
−12.566 ± j11107, −12.566 ± j11107
fig.3.9(e) −37.762 ± j7639.4, −37.762 ± j7639.4, −24.64 ± j3298.6, −24.64 ± j3298.6,
−15.188 ± j16584, −15.188 ± j16584
Table 3.2: Eigenvalues of the systems given in fig. 3.9

Chapter 3 31
Distributions of the Characteristic Eigenvalues of the Study Network I and Its Vari-
ants
The CE of the SNI are plotted in fig. 3.11(a). It is worth noting that all the CE of the
network in fig. 3.6 have similar real parts though the imaginary parts vary considerably.
Further, these real part values match those given in fig. 3.8(d). In order to get some more
general conclusions about the distribution of the CE, more cases are studied, which are
variants of the system in fig. 3.6. The variations include both the network topology and
the line parameters.
The variant networks are shown in fig. 3.10, in which, the line between nodes 824 and
828 is removed, and a new line is added between nodes 826 and 834 in fig. 3.10(a) and
between nodes 814 and 858 in fig. 3.10(b) respectively.
Figure 3.10: Variants of the SNI
The two topologies in fig. 3.10 will be referred to as ’T2’ and ’T3’, while the topology
in fig. 3.6 will be referred to as ’T1’. For line parameters, the variations as well as the
32 Chapter 3
original one can also be separated as three groups:
• For the first group (later referred to as ’P1’), the line parameters are given in ap-
pendix B. If the length of one of the lines is 1km, its CE is −164.8 + j166850, as
given in fig. 3.8.
• For the second one (later referred to as ’P2’), the parameters of the lines 808-812,
816-824 and 860-836 are altered as: Rn = 2 × R, Ln = L/2, Cn = C, where R, L, C
have the same values given in appendix B, while Rn , Ln , Cn the altered parameters
of the lines. If the length of one of the altered lines is 1km, its CE is −659.22 +
j235960.
• For the third group (later referred to as ’P3’), besides the alteration of the lines
808-812, 816-824 and 860-836 as in the second group, lines 830-854, 854-856 and
854-852 are also changed as: Rn = R/2, Ln = 2 × L, Cn = C. If the length of one of
the altered lines is 1km, its CE is −41.2 + j117980.
Detailed information of each studied case is summarized in table 3.3.
Case No. Network Topology Line Parameters Plot of Characteristic Eigenvalues

1 T1 P1 fig.3.11(a)
2 T2 P1 fig.3.11(b)
3 T3 P1 fig.3.11(c)
4 T1 P2 fig.3.11(d)
5 T2 P2 fig.3.11(e)
6 T3 P2 fig.3.11( f )
7 T1 P3 fig.3.11(g)
8 T2 P3 fig.3.11(h)
9 T3 P3 fig.3.11(i)
Table 3.3: Detailed information of the SNI and the variants
From fig. 3.11(a)-(c), it can be seen that though the network topology changes, the real
parts of the CE remain unchanged, which match the one given in fig. 3.8(d). This indi-
cates that, compared to network topology, line parameters predominantly determine the
Chapter 3 33
real parts of the CE. This observation can be verified by investigating fig. 3.11(d)-(i). Us-
ing fig. 3.11(d)-( f ) as an example, in which the distribution of the CE are similar though
the differences among them are distinct; while the comparison of fig. 3.11(a)-(c) to fig.
3.11(d)-( f ) seems to reveal that the CE are subjected to a large change. Hence, though
the topology has an influence on the CE, the parameters of the line are of much greater
influence.
A more instructive principle about the relation between the line parameters and the CE of
the network can be formed as: assuming a network consists of M different types of lines,
then the maximum and minimum values of the real parts of their CE span a range, which
covers the real part values of the CE of the network. An explanation is given based on
case 9. In it, the network has 3 different types of lines. The real parts of their CE are
-41.2, -164.8 and -659.22 respectively. Their maximum and minimum values of the real
parts span a range from -41.2 to -659.22, which totally covers the range spanned by the
real parts of the CE of the given network, which from fig. 3.11(i) can be noted as from
-44.9 to -526.2.
34 Chapter 3
Figure 3.11: Distributions of the characteristic eigenvalues of the SNI and the variants
Chapter 3 35
Distributions of the Characteristic Eigenvalues of the Study Network II
The CE of the SNII are plotted in fig. 3.12. Meanwhile, for ease of comparison, the CE
of the cables used in the SNII are given in table 3.4.
Figure 3.12: Distributions of the characteristic eigenvalues of the SNII
Type of cable CE Type of cable CE

NA2XS(F)2Y 3×1×185 -160.01 + j138820 NKBA 3×70 -409.31 + j191350
NAKBA 3×240 -229.44 + j161700 NAKBA 3×95 -507.37 + j184290
NAKBA 3×150 -350.92 + j172360 NAKBA 3×70 -673.63 + j191350
Table 3.4: Characteristic eigenvalues of the cables used in the SNII†
The results indicate that as for the cables of the SNII, the maximum and minimum values
of the real parts of their CE span a range from -160.01 to -673.63, which totally covers
the range spanned by the real parts of the CE of the SNII, which from fig. 3.12 can be
noted as from -160.01 to -513.27. This conforms to the conclusions made in the previous
part.
† When calculating the CE, all the cables are assumed to have a length of 1km.
36 Chapter 3
3.4.2 Analysis of the Characteristics of Distribution Network Dy-

namics
Neplan [46] offers comprehensive information about the available lines/cables for the dis-
tribution systems with a voltage level from 1kV to 20kV in Europe. In its databank, the
parameters of 361 types of lines/cables are given. These data are fully studied and the
results reveal that, the maximum real part value of the CE of these lines/cables is -69.3,
while the minimum value of that is -33055. Thus it can be expected that, the eigenvalues
of a real distribution network in Europe will have real parts sited within the range from
-69.3 to -33055. This means that a dynamic procedure caused by a certain disturbance
will die down to 5% of its initial value within 0.0001s to 0.0433s. Though without proof
of enough statistics, it is sensible to judge that normally the time intervals between two
disturbances will be much bigger than 0.0433s. A bit work of clarifying the meaning of
’disturbance’ used here may be helpful. In this thesis, the stochastic properties of DG are
also deemed as kind of disturbances, though practically the alteration of its status may be
the result of some normal operations, for example, the switching on/off of it.
From aforesaid analysis it can be seen that though apparently the disturbances are stochas-
tic series, since the corresponding dynamic procedures last only quite a short period, the
impact of each individual of the series to the electrical network is nearly independent.
Besides, under this condition, it is still proper to adopt the steady state or quasi steady
state assumptions to model the distribution network with DG.
3.4.3 Conclusions and Comments
This chapter focuses on studying the attributes of the electrical network dynamics of dis-
tribution systems. By investigating the eigenvalues of the state space model of the elec-
trical network, the following conclusions are drawn:
• There is a close relation between the eigenvalues of the network under abc and dq0
reference frame. Assuming a symmetric network has 3N state variables, thus it
has 3N eigenvalues. Among them, N pairs of the complex-conjugate/double-real
eigenvalues under abc frame, relate only to the network dq components under dq0
Chapter 3 37
frame. These eigenvalues under dq0 frame can be gotten by shifting the counter-
parts under abc frame with frequency of ±ω0 . The remaining N eigenvalues will
keep unchanged under both of the two frames.
• The parameters of the lines/cables of the electrical network are the predominant
factors to determine the eigenvalues of the network. Although the network topology
has an influence on the eigenvalues, the parameters of the lines/cables are of much
greater influence.
• The distribution of the eigenvalues of the network is bounded by the relevant values
of its lines/cables. Assuming a distribution network consists of M different types
of lines, then the maximum and minimum values of the real parts of their CE span
a range, which covers the total interval spanned by the real parts of the CE of the
network.
• In a real distribution network in Europe, the eigenvalues will have the real parts
within the range from -69.3 to -33055. This is based on a thorough investigation of
the data of the lines/cables given in Neplan [46].
• The impacts of the apparently successive stochastic disturbances to the network

can be studied independently and the nodal admittance matrix is acceptable as the
model of a distribution network. The investigation on the eigenvalues of a distribu-
tion network reveals that the electrical dynamics die off very fast. Hence a steady
state assumption for modeling the distribution network is acceptable. The appar-
ently successive stochastic disturbances have nearly independent influences on the
network.
So far, some instructive conclusions are obtained. However, some limitations about the
results exist as well. They are:
• The eigenvalues of set λabcN are not included in the discussion of section 3.4.1.
However, it can be reasonably expected that the concerned conclusion is also appli-
cable to λabcN . That is, assuming the network consists of M different types of lines,
then the maximum and minimum values of the real parts of their CE in λabcN span
38 Chapter 3
network in λabcN .
• Frequency dependence of line parameters is not considered when calculating the

eigenvalues. Thus the results may not match the practical observations. A possibly
improved approach is: 1. assuming the system is stepwise frequency independent;
2. forming the state matrix within this frequency band ([47]-[51] are some refer-
ences for modeling frequency dependent lines); 3. calculating the eigenvalues, and
those whose imaginary parts conform to the predetermined frequency band will
belong to the set of the eigenvalues of the real system; 4. step 1 to 3 can be exe-
cuted iteratively until all the eigenvalues have been found. During this procedure,
it is important to note that as frequency keeps rising, the π equivalent model for
lines/cables will not be proper even for distribution networks, in which the length of
the line/cable normally is short. Under this situation, distributed parameter models
can be adopted instead. However, the network dynamic order and thus the quantity
of the eigenvalues will increase.
• All the examples and studied systems in this chapter are based on radial network
topology. Thus it is open for further study if the introduced results are also applica-
ble to meshed networks. However, a positive answer can be reasonably expected.
Chapter 4
Impacts of DG on the Small Signal

Stability of Distribution Systems
Small signal stability concerns the ability of a system to maintain stability under small
disturbances. Modal analysis, as a common tool to investigate the small signal stability,
is adopted in this chapter to investigate the influence of DG on the small signal stability
of distribution systems.
In this chapter, section 4.1 gives a brief view of the implementation of the small signal
stability analysis in power systems. Section 4.2 presents the models for all the relevant
system components. Section 4.3 contains simulations, results and comments.
The work of this chapter shows that the integration of DG will not lead to any small
signal instability problem, even if the penetration level of DG is high. Besides, proper
shunt reactive compensation, such as the use of STATCOM, improves the small signal
stability of the system as well as the penetration level of DG.
40 Chapter 4
4.1 Implementation of the Small Signal Stability Analy-

sis in Power Systems
4.1.1 Basic Concepts of Modal Analysis
For a power system, if the dynamics of the electrical network are also included into the
small signal stability analysis, the whole system will be described by an ODE model.
However, as discussed in chapter 3, since the dynamics of the electrical network die out
very fast, these dynamics can be neglected which hence leads to the assumption that the
electrical network is always in a steady state. Under this assumption, a DEA model of the
whole system will be built instead of an ODE one.
Generally, a DAE model of a power system has the form as
ẋ = f (x, v, u)
0 = g(x, v) (4.1)
where
x is a vector of state variables;
v a vector of algebraic variables;
u a vector of disturbances or control inputs;
f represents the nonlinear dynamic characteristics of system components;
g the nonlinear relations between x and v.
A linearized model around an operating point can be deduced from eqn. (4.1) as
∆ẋ = J f x ∆x + J f v ∆v + J f u ∆u
0 = Jgx ∆x + Jgv ∆v (4.2)
where J f x , J f v , J f u , Jgx and Jgv are the respective Jacobian matrices.
Eliminating the algebraic variables v in eqn. (4.2) leads to
∆ẋ = A∆x + B∆u (4.3)

Chapter 4 41
−1 J , B = J .
where A = J f x − J f v Jgv gx fu
In eqn. (4.3), A is the state matrix. The small signal stability can be evaluated by in-
vestigating its eigenvalues.
Mode shape and eigenvectors [21]

Assuming A a n × n matrix, thus it has n eigenvalues:
λ = λ1 , λ2 , ...λn
For each eigenvalue λi , its associated n-column right eigenvector Φi satisfies the following
equation:
AΦi = λi Φi (4.4)
Similarly, the so-called left eigenvector Ψi contains n rows and satisfies the equation:
Ψi A = λi Ψi (4.5)
Physically, the right eigenvector describes how each mode of oscillation is distributed
among the systems states and it is called mode shape. The left eigenvector determines the
combination of the state variables on a specific mode.
Participation matrix [21]

The participation matrix P combines the right and left eigenvectors to measure the relation
between the state variables and the modes.
P = [p1 p2 ... pn ] (4.6)
with    
p1i φ1i ψi1
   
 p   φ ψ 
 2i   2i i2
pi =  =


 ...   ... 
   
pni φni ψin
where
φki is the element on the kth row and ith column of the modal matrix Φ, with Φ = [Φ1 Φ2
... Φn ];
42 Chapter 4
ψik is the element on the ith row and kth column of the modal matrix Ψ, with Ψ = [ΨT1 ΨT2
... ΨTn ]T ;
pki = φki ψik is the participation factor measuring the relative participation of the kth state
variable in the ith mode, and vice versa.
Damping ratio [21]

The damping ratio ξ determines the decaying rate of the amplitude of an oscillation. For
a complex-conjugate pair of eigenvalues
λ = σ ± jω (4.7)
the corresponding damping ratio is defined as

−σ
ξ= √ (4.8)
σ2 + ω2
4.1.2 Selection of a Common Reference Frame

In a multimachine system, it is normal to refer all the models to a CRF. The determination
of the CRF is presented in section 3.1.1.
4.1.3 Some Specifics of the Power System Model

The formulation of the state matrix of a power system involves steps of linearizing com-
ponent models and eliminating all the algebraic variables. Various components and their
large quantities make this procedure very complex. In order to form the state matrix A in
a systematic way, in this thesis the algebraic variables in eqn. (4.1) are specifically cho-
sen as the amplitude and argument of the nodal voltages, while the component models are
coupled with the network equations via the principle of power conservation. This means
that in eqn. (4.1), the nonlinear equations g can be rewritten as
0 = g1 (x, v) + g2 (v) (4.9)
where
g1 describes the quantities of the power generated/consumed by system components, such
as generators, loads, FACTS devices, etc;
g2 is the power transferred by the network, actually, which has the same form as in a
Chapter 4 43
traditional power flow calculation.
The principle of the power conservation is reflected by the sum of g1 and g2 in eqn.
(4.9) equal to zero.
4.2 Models of Power System Components
4.2.1 Equivalent of the Transmission System

In this thesis, since more attention is paid to the influence of DG on the distribution sys-
tem, the transmission system is modeled by a synchronous generator with a large capacity.
The classical model is adopted here for the reason that the impact of the operation status
of the distribution system to the transmission system is not strong hence the equivalent
synchronous generator need not to be modeled in detail. The corresponding mathematical
model is
δ˙ = ω0 (ω − 1)
ω̇ = (Pm − Pe − D(ω − 1))/M (4.10)
Where Pm is the mechanical power, Pe the electrical power, D the damping coefficient and
M the inertia constant.
The electrical power Pe can be calculated as
Pe = (vq + ra iq )iq + (vd + ra id )id (4.11)
Meanwhile, the relations between the voltages and the currents are
0 = vq + ra iq − e0q + xd0 id
0 = vd + ra id − xd0 iq (4.12)
where e0q is the q-axis transient voltage, ra the stator resistance, and xd0 the d-axis transient
reactance.
In eqn. (4.11) and (4.12), vd and vq are voltage components referred to the individual
44 Chapter 4
dq0 frame of the synchronous generator. They are connected with the stator voltage un-
der the CRF via
vd = V sin(δ − θV )
vq = V cos(δ − θV ) (4.13)
where V is the amplitude of the stator voltage, θV the angle referred to the CRF.
Similarly, id and iq are current components referred to the individual dq0 frame of the
synchronous generator. However, it is not necessary to formulate some equations similar
to eqn. (4.13), to relate id and iq to the stator current under the CRF, since they have
already been connected to the algebraic variables vd and vq through eqn. (4.12).
The power injections from the generator are
P = vd id + vq iq
Q = vq id − vd iq (4.14)
4.2.2 Network
As analyzed in chapter 3, since the dynamics of a distribution network die out very fast,
when studying the small signal stability, it is proper to assume that the network is al-
ways in a steady state, or more accurately, a quasi steady state. Thus the corresponding
mathematic model of the network is a nodal admittance matrix.
4.2.3 Constant Speed Wind Turbine Generation System

The DG in this chapter are specifically designated as constant speed wind turbine genera-
tion systems. The reason is that: compared to other DG technologies, the mathematically
mature models of the wind generation systems (WGS) are available; further, among all
the available WGS, the constant speed wind turbine generation system has the biggest
interaction with the connected power system.
The model of the constant speed wind turbine generation system used in this work is
presented in [31], which stems from the original one proposed in [4].
Chapter 4 45
The Mechanical Power Extracted by Wind Turbine from the Wind
The mechanical power extracted by wind turbine from the wind can be calculated as
ρ
Pw = c p (λ)Ar v3w (4.15)
2
where, ρ is the air density, c p the performance coefficient or power coefficient, Ar the area
covered by the wind turbine rotor and λ the tip speed ratio (the ratio between the blade tip
speed vt (m/s) and the wind speed vw (m/s) at the hub height upstream the rotor):
vt 2Rωwr
λ= = ηGB (4.16)
vw pvw
with ηGB the gear box ratio, p the number of poles of the induction generator, R the rotor
radius, and ωwr the mechanical angular speed of the wind turbine.
Additionally, the following equation is used to approximate the c p curve:
125 − 16.5
c p = 0.44( − 6.94)e λi (4.17)
λi
1
with λi = 1
+0.002
.
λ
Shaft Model
A two mass representation has been adopted to describe the dynamic behavior of the
mechanical driving subsystem which delivers the mechanical power from the rotor to the
induction generator:
Twr − Ks γ
ω̇wr =
2Hwr
Ks γ − Te
ω̇m =
2Hm
γ̇ = 2π f (ωwr − ωm ) (4.18)
where, f is the nominal grid frequency, γ the angular displacement between the two ends
of the shaft, Ks the shaft stiffness, Twr the torque, Hwr the inertia constant of the wind tur-
bine, Te and Hm the corresponding values of the generator, and ωm the mechanical angular
speed of the generator.
46 Chapter 4
In equation (4.18), the mechanical torque Twr and the electrical torque Te are given by:
Pw
Twr =
ωwr
Te = e0r ir + e0m im (4.19)
while the explanation of e0r , ir , e0m and im is in the following paragraph.
Generator Model
The generator used in a constant speed generating system is a squirrel cage induction
generator. Single cage representation for the rotor is adopted. When the stator transients
are neglected, the model of the generator has the order of 2 and can be expressed as
follows:
e0r − (x0 − x0 )im

ė0r = ω0 (1 − ωm )e0m −
T00
e0 + (x0 − x0 )ir
ė0m = −ω0 (1 − ωm )e0r − m (4.20)
T00
with algebraic constraints as
e0r − vr = rs ir − x0 im
e0m − vm = rs im + x0 ir (4.21)
where x0 , x0 and T00 can be obtained from the generator parameters of rotator resistance
rR , rotator reactance xR , stator reactance xS and magnetizing reactance xm as:
x0 = xS + xm
xR xm
x0 = xS +
xR + xm
x R xm
T00 = (4.22)
ω0 rR
and vr , vm , ir and im are the projections of the stator voltage and current to the CRF.
The following relations between vr , vm and the stator voltage hold
vr = V sin(−θV )
vm = V cos θV (4.23)
where V and θV are the amplitude and argument of the stator voltage in the CRF.
Chapter 4 47
Electric Power Generated by the Wind Energy Conversion
The generated power is given by
P = vr ir + vm im
Q = vm ir − vr im + bc (v2r + v2m ) (4.24)
where bc is the compensating capacitor conductance.
4.2.4 STATCOM
Based on the works in [52]-[57] and [7], an improved STATCOM model is developed and
presented in this section.
Figure 4.1: Electrical diagram of STATCOM
Pulse width modulation (PWM) is assumed to control the STATCOM via two principles:
(1) by changing the ratio m to adjust the magnitude of the ac voltage and thus the exchange
of the reactive power with the network; (2) by changing the phase angle ϕ between the
network voltage and the voltage generated by the converter to feed the appropriate power
from the network to balance the power loss of the STATCOM.
Before developing the STATCOM model, it is important to note that the per unit system
adopted here is:
48 Chapter 4
• for DC variables and AC variables in abc reference frame, the base voltage is the
rated line-to-line voltage of the STATCOM.
• for variables in dq0 reference frame, the base voltage is the peak value of the rated
line-to-neutral voltage of the STATCOM.
In both cases, the base power is the rated capacity of the STATCOM. Thus the following
relations are valid:
√
Sb = 3Vabc Iabc
√
Zabc = Vabc / 3Iabc
Cabc = 1/ω0 Zabc (4.25)
and
3
Sb = Vdq0 Idq0
2
r
3
Vabc = Vdq0 (4.26)
2
where, all the given values are the base quantities in the respective reference frames.
In fig. 4.1, the DC side dynamics satisfy

Idc
V̇dc = (4.27)
Cdc
Inserting eqn. (4.25) into eqn. (4.27) leads to the per unit model
ω0 I¯dc
V̄˙dc = √ (4.28)
3C̄dc
where the superbar denotes the per unit variables.
Since the instant power at the AC and DC side must be equal, this results in
2
Vdc
Vdc Idc + = Vaca Iaca +Vacb Iacb +Vacc Iacc (4.29)
Rdc
where Vaca and Iaca are the AC voltage and AC current of phase a respectively, and so on.
Transforming eqn. (4.29) to dq0 frame and inserting eqn. (4.26) into it leads to
1 V̄ 2
√ V̄dc I¯dc + dc = V̄acd I¯acd + V̄acq I¯acq (4.30)
3 R̄dc
Chapter 4 49
The AC voltage at the converter is
V̄ac = mkV̄dc 6 ϕ (4.31)
where m is the modulation ratio defined by PWM, k a ratio depending on the converter
q
structure [58], in this work k = 38 .
It is noted in advance that in the following equations, all the variables are expressed in per
unit quantities. However, the superbar is omitted to make the equations concise.
Referring eqn. (4.31) to the CRF results in
Vacd = mkVdc sin(−ϕ)

Vacq = mkVdc cos ϕ (4.32)
Inserting eqn. (4.32) into eqn. (4.30) leads to

√
√ 3Vdc
Idc = 3mk(Iacd sin(−ϕ) + Iacq cos ϕ) − (4.33)
Rdc
Substituting eqn. (4.33) into eqn. (4.28) leads to
ω0 ω0Vdc
V̇dc = mk(Iacd sin(−ϕ) + Iacq cos ϕ) − (4.34)
Cdc RdcCdc
Assuming ~Iac = Iac 6 θI (the superarrow denotes a phasor), this means
~Iac = Iac cos θI + jIac sin θI (4.35)
Referring ~Iac to the CRF leads to
Iacd = Iac sin(−θI )

Iacq = Iac cos θI (4.36)
Combining eqn. (4.35) with (4.36) results in
Iacd = −imag(~Iac )
Iacq = real(~Iac ) (4.37)
50 Chapter 4
Meanwhile, assuming ~Vt = Vt 6 θ leads to

~ ~
~Iac = Vt − Vac
Zt
= (Vt cos θ −Vac cos ϕ)Gt − (Vt sin θ −Vac sin ϕ)Bt
+ j((Vt cos θ −Vac cos ϕ)Bt + (Vt sin θ −Vac sin ϕ)Gt ) (4.38)
Combining eqn. (4.37) with (4.38) leads to
Iacd = −(Vt cos θ −Vac cos ϕ)Bt − (Vt sin θ −Vac sin ϕ)Gt
Iacq = (Vt cos θ −Vac cos ϕ)Gt − (Vt sin θ −Vac sin ϕ)Bt (4.39)
Inserting (4.39) into (4.34) results in

ω0 ω0Vdc
V̇dc = mk(Vt sin(ϕ − θ)Bt + (Vt cos(ϕ − θ) − mkVdc )Gt ) − (4.40)
Cdc RdcCdc
In this work, two PI controllers shown in fig. 4.2 are adopted to control the m and ϕ in
eqn. (4.31).
Figure 4.2: Controllers of STATCOM
It is worth noting that in fig. 4.2(b) the signs of Vdcre f and Vdc are just opposite to those
given in [53] and [54]. This can be explained with the help of fig. 4.3. Normally, Vac lags
behind Vt since the STATCOM will absorb active power from the network. This active
power is used partly to compensate the STATCOM power losses and partly to charge the
capacitor to increase its voltage. It is assumed that initially the STATCOM is under the
steady state, thus Vdc equals Vdcre f . Due to a disturbance, Vdc increases a bit. The dif-
ference between Vdc and Vdcre f will generate a positive signal ∆ϕ which makes Vac rotate
0
counterclockwise a bit (from Vac to Vac ), thus decreases the angle gap between Vt and Vac .
This will reduce the active power fed by the network to the STATCOM, and consequently
Chapter 4 51
decrease the voltage Vdc . Similarly, when a disturbance makes Vdc decrease a bit, analysis
shows that the controller will try to make an increment of it. This indicates that the design
in fig. 4.2(b) is stable.
Figure 4.3: Explanation of the structure of the DC voltage controller in fig. 4.2
In the following the mathematical models of the STATCOM controllers are developed.
The AC voltage controller in fig. 4.2(a) is defined by
Kacp s + Kaci
m= (Vtre f −Vt ) (4.41)
s
Introducing
1
x1 = (Vtre f −Vt )
s
leads to the time domain representation
ẋ1 = Vtre f −Vt

m = Kacp (Vtre f −Vt ) + Kaci x1 (4.42)
The DC voltage controller in fig. 4.2(b) is given by
Kdcp s + Kdci
ϕ= (−Vdcre f +Vdc ) (4.43)
s
Introducing
1
x2 = (−Vdcre f +Vdc )
s
52 Chapter 4
leads to
ẋ2 = −Vdcre f +Vdc

ϕ = Kdcp (−Vdcre f +Vdc ) + Kdci x2 (4.44)
The power absorbed by the STATCOM from the network is given by
P = Vt2 Gt − mkVt Vdc Gt cos(θ − ϕ) − mkVt Vdc Bt sin(θ − ϕ)

Q = −Vt2 Bt + mkVt Vdc Bt cos(θ − ϕ) − mkVt Vdc Gt sin(θ − ϕ) (4.45)
Eqn. (4.40), (4.42), (4.44) and (4.45) constitute the complete model of the STATCOM.
In order to include the developed model of the STATCOM into the modal analysis, a
linearization of the model is necessary.
The linearized model can be written as
∆ẋ = J f x ∆x + J f v ∆v + J f u ∆u
∆g1 = Jg1 x ∆x + Jg1 v ∆v
where x is the state variable vector, which includes [Vdc , x1 , x2 ]T

v the algebraic variable vector, v = [θ, Vt ]T
u the input variable vector.
The Jacobian matrices are defined as
ω0 2 2 ω0 ω0 ω0
J f x (1, 1) = − m k Gt − +( mkVt cos(ϕ − θ)Bt − mkVt sin(ϕ − θ)Gt ) ∗ Kdcp
Cdc RdcCdc Cdc Cdc
ω0 ω0 ω0 2
J f x (1, 2) = ( kVt sin(ϕ − θ)Bt + kVt cos(ϕ − θ)Gt − 2 mk Vdc Gt ) ∗ Kaci
Cdc Cdc Cdc
ω0 ω0
J f x (1, 3) = ( mkVt cos(ϕ − θ)Bt − mkVt sin(ϕ − θ)Gt ) ∗ Kdci
Cdc Cdc
J f x (3, 1) = 1
Chapter 4 53
ω0 ω0
J f v (1, 1) = − mkVt cos(ϕ − θ)Bt + mkVt sin(ϕ − θ)Gt
Cdc Cdc
ω0 ω0
J f v (1, 2) = mk sin(ϕ − θ)Bt + mk cos(ϕ − θ)Gt +
Cdc Cdc
ω0 ω0 ω0 2
−( kVt sin(ϕ − θ)Bt + kVt cos(ϕ − θ)Gt − 2 mk Vdc Gt ) ∗ Kacp
Cdc Cdc Cdc
J f v (2, 2) = −1
Jg1 x (1, 1) = −mkVt Gt cos(θ − ϕ) − mkVt Bt sin(θ − ϕ) +
(−mkVt Vdc Gt sin(θ − ϕ) + mkVt Vdc Bt cos(θ − ϕ)) ∗ Kdcp
Jg1 x (1, 2) = (−kVt Vdc Gt cos(θ − ϕ) − kVt Vdc Bt sin(θ − ϕ)) ∗ Kaci
Jg1 x (1, 3) = (−mkVt Vdc Gt sin(θ − ϕ) + mkVt Vdc Bt cos(θ − ϕ)) ∗ Kdci
Jg1 x (2, 1) = mkVt Bt cos(θ − ϕ) − mkVt Gt sin(θ − ϕ) +
(mkVt Vdc Bt sin(θ − ϕ) + mkVt Vdc Gt cos(θ − ϕ)) ∗ Kdcp
Jg1 x (2, 2) = (kVt Vdc Bt cos(θ − ϕ) − kVt Vdc Gt sin(θ − ϕ)) ∗ Kaci
Jg1 x (2, 3) = (mkVt Vdc Bt sin(θ − ϕ) + mkVt Vdc Gt cos(θ − ϕ)) ∗ Kdci
Jg1 v (1, 1) = mkVt Vdc Gt sin(θ − ϕ) − mkVt Vdc Bt cos(θ − ϕ)
Jg1 v (1, 2) = 2Vt Gt − mkVdc Gt cos(θ − ϕ) − mkVdc Bt sin(θ − ϕ) +
(kVt Vdc Gt cos(θ − ϕ) + kVt Vdc Bt sin(θ − ϕ)) ∗ Kacp
Jg1 v (2, 1) = −mkVt Vdc Bt sin(θ − ϕ) − mkVt Vdc Gt cos(θ − ϕ)
Jg1 v (2, 2) = −2Vt Bt + mkVdc Bt cos(θ − ϕ) − mkVdc Gt sin(θ − ϕ) +
−(kVt Vdc Bt cos(θ − ϕ) − kVt Vdc Gt sin(θ − ϕ)) ∗ Kacp
4.2.5 Loads
The static load model is adopted here as to represent the active power as constant current
component and the reactive power as constant impedance component.
4.3.1 Impacts of WGS on the Small Signal Stability of the SNI
When investigating the impact of WGS on the small signal stability of the SNI, the fol-
lowing factors are expected to have a relatively large influence on the results:
54 Chapter 4
• Load Level —- A high load level is selected as corresponding to the load values in
appendix B; while a low load level corresponds to 20% of the demand of the high
load level.
• Characteristics of the Connection of the Distribution System to the Transmission

System —- A strong connection is selected as corresponding to the line segment
at 800-802 in fig. 3.6 with the specified length as in appendix B; while a week
connection is obtained by stretching this line 10 times longer to 7.864 km.
• WGS penetration level —- Though there are no practical evidences, it can be rea-
sonable to expect that the more WGS connected to the network, the higher their
influences on the system operation. Thus it is assumed the WGS achieves its maxi-
mum allowable penetration level.
• WGS generation level —- The generation level is mainly determined by the wind
speed, which in this work varies discretely from 0.4 p.u. to 1.0 p.u. rated wind
speed (approximately corresponding to no generation to full generation) with a step
of 0.1 p.u..
All the scenarios are classified into 4 groups according to the combinations of the above
mentioned factors. Details are given in table 4.1.
Group 1 Group 2 Group 3 Group 4

load level high low high low
connection characteristic of
strong weak
distribution to transmission systems
WGS penetration level maximum
WGS generation level stepwise changing from no to full generation
Table 4.1: Detailed information of the studied scenarios
Determination of the Maximum Allowable Penetration Level of WGS
In this work, the maximum allowable penetration level of WGS is evaluated by the max-
imum number of WGS allowed to connect to the distribution network, which means that
Chapter 4 55
for a given distribution system, if the maximum allowable number of WGS is n, then for
all the possible scenarios, the system must be in a normal operation state; however, if
there are totally n + 1 WGS integrated to the network, then for some scenarios the system
will be subjected to abnormal operation problems, which can be a steady state voltage
violation, or a small signal instability problem, etc.
The penetration level is preliminarily determined by considering the fact that the steady
state voltage of the network must stay within a fixed predetermined range. Here the al-
lowable steady state voltage variation is adopted equal to ±10% (0.9/1.1 p.u.). Later, the
feasibility of the obtained penetration level will be verified by checking the status of the
small signal stability.
In the following parts, the simulation procedures and results for the scenarios of group
1 and 2 are given in detail. At the end of this section, a brief view of the results for the
scenarios of group 3 and 4 are available.
The maximum allowable numbers of the WGS connected to each node for group 1 and
2 are determined with the help of power flow computations. The results reveal that, the
higher the generation level of the WGS, the lower the allowable penetration level of the
WGS.
Fig. 4.4(a) and 4.4(b) present the results for the scenarios of group 1 and 2 respectively.
To make the figures clear, only the results with the assumption that all the WGS are at the
full power generation are shown.
From fig. 4.4 it can be seen when the WGS is connected to different nodes, the maximum
allowable numbers of it vary quite a lot. It is worth emphasizing again, that these results
are obtained according to the criterion of the nodal steady state voltage.
As stated in table 4.1, when investigating the small signal stability, the integration of
the WGS is assumed to achieve the maximum level. In order to determine the specific
number of the integrated WGS by referring to the results given in fig. 4.4, the following
56 Chapter 4
Figure 4.4: The maximum allowable numbers of the WGS connected to each node for
group 1 and 2
Chapter 4 57
considerations are relevant:
• In this thesis, since the objective is to evaluate the impacts of the DG rather than
only the WGS on system operations in a more general sense, hence the WGS is
assumed to be allowed to connect to the network arbitrarily, which may give a more
suitable approximation to the potential operation of the DG in the future thus lead to
a more general conclusion. This assumption implicitly puts forward a requirement
as, to determine the number of the WGS from the presently available results given
in fig. 4.4, it must be note that under all the possible distributions of the WGS, the
steady state voltage must be in the fixed range. This requirement can be met by
choosing the number of the WGS as the minimum one given in the figure, which
for fig. 4.4(a) is 11, for fig. 4.4(b) is 7.
• Scenarios of group 1 and 2 should have the same maximum allowable penetration
level of WGS, since they are based on a physically identical network.
The above considerations result in the conclusion that for the scenarios of group 1 and 2,
the maximum number of WGS allowed to connect to the network is 7.
Impacts of WGS on the Small Signal Stability
As above mentioned, for the scenarios of group 1 and 2 the studied network allows at
most 7 WGS to be connected. Before moving forward to deeper analyses, a general view
about the possible distribution of the eigenvalues of the studied system may be helpful.
Table 4.2 gives a sample of the eigenvalues of the system with conditions as: the 7 WGS
are connected to node 838; the wind blowing on them has speed as 0.4, 0.5, 0.6, 0.7, 0.8,
0.9 and 1.0 respectively.
A sensible model of wind speed for different WGS will consider the correlation among
them. This is necessary when the long term system operations with WGS are investi-
gated. However, this work focuses on studying the small signal stability at some certain
moments, at which the system status is got by snapshooting the practical system state.
Under this condition, the simplification as assuming the wind speeds at different spots are
58 Chapter 4
independent, as adopted in the above discussion, is acceptable.
Damping Damping
No. Eigenvalue Ratio (%) No. Eigenvalue Ratio (%)
1 -17.222 100 17, 18 -8.415 ± j43.412 19.03
2 -17.342 100 19, 20 -8.438 ± j43.405 19.083
3 -17.478 100 21, 22 -3.918 ± j33.325 11.677
4 -17.604 100 23, 24 -0.462 ± j4.232 10.845
5 -17.746 100 25, 26 -0.458 ± j4.232 10.767
6 -17.696 100 27, 28 -0.449 ± j4.228 10.563
7 -10.771 100 29, 30 -0.428 ± j4.224 10.074
8 -0.789 100 31, 32 -0.395 ± j4.218 9.327
9, 10 -8.653 ± j43.129 19.67 33, 34 -0.354 ± j4.209 8.387
11, 12 -8.600 ± j43.205 19.523 35, 36 -0.311 ± j4.196 7.398
13, 14 -8.540 ± j43.29 19.354 37 0 -
15, 16 -8.482 ± j43.362 19.198
Table 4.2: A sample of the eigenvalues for a specific scenario
As shown in table 4.2, when 7 WGS are connected to the studied network, the whole sys-
tem totally has 37 eigenvalues. In them, the zero eigenvalue numbered 37 appears just due
to the lack of uniqueness of absolute rotor angle [21], thus requires no further attention.
However, it can be seen that the pair of eigenvalues numbered 35 and 36, have the min-
imum absolute values of the real part as well as the minimum damping ratio. Thus they
may play a dominant role in the dynamic responses of the system. Therefore the impacts
of WGS on the small signal stability can be investigated by studying the evolution of this
pair of eigenvalues when other conditions, such as the locations, the generation levels of
the WGS are changed. For ease of later reference, this pair of eigenvalues is named as
’Critical Eigenvalues’ (CRE).
To investigate the influences of the integrated WGS on the small signal stability, factors
as the locations and the generation levels of the WGS are expected to be most relevant.
The total number of the possible scenarios makes an exhaustive study quite hard if not
Chapter 4 59
impossible. Thus another analysis procedure consisting of three steps is proposed as fol-
lows:
• Firstly, evolutionary algorithm is used to form a conjecture about the relation be-
tween the concerned factors and the small signal stability.
• Then, the obtained conjecture is exploited as a direction for seeking the potentially
worst scenario.
• Finally, lots of randomly generated scenarios are analyzed to verify the results.
When the correctness of the results is justified, some general conclusions are drawn
from them.
The details of the procedure are:
Step 1, applying evolutionary algorithm to identify the causality between the concerned
factors and the potential problems of the small signal instability.
Evolutionary algorithm [59]-[64] has been successfully used for power system optimiza-
tions in many aspects. Since the interest of this work is to identify if a high penetration
level of DG will induce the problems of the small signal instability, to find the potential
scenario of which the CRE have the minimum damping ratio (compared to the damping
ratio of the CRE of other scenarios), thus will be meaningful and important. The seeking
of such a scenario can be proposed as a optimization problem, i.e., to find a combination
of the concerned factors to achieve the minimum damping ratio of the CRE. Thus evolu-
tionary algorithm is applicable in this case.
The structure of the evolutionary algorithm used in this work is shown in fig. 4.5.
The details of the implementation of the algorithm are as follows:
• In the stage of initialization, the population is set to have 50 individuals. Since

there are 7 WGS under the consideration, further, for each WGS the location and
the generation level are variables to be determined, thus each individual is coded to
include 14 units, in which, the first 7 successive units indicate the locations of the
60 Chapter 4
Figure 4.5: Structure of the evolutionary algorithm [65]
WGS, while the last 7 successive units record the wind speed. A schematic of such
an individual is given in fig. 4.6.
Figure 4.6: Schematic of a individual
• The objective function is to get the minimum value of the damping ratio of the CRE
under all calculated scenarios.
• The optimization criteria, which serves actually as a convergence condition, is

specifically chosen as the maximum generation number under which the value of
the objective function keeps unchanged is more than 100.
• The fitness for the stage of selection is calculated through the following function
1
fi = ( )2 (4.46)
10 × DRi
where fi is the fitness of the ith individual, DRi the damping ratio of the CRE of the
ith individual.
Besides, roulette-wheel selection is adopted as the selection scheme.
• The possibility of recombination is 0.5.

Chapter 4 61
• The possibility of mutation is 0.2.
Fig.4.7(a) and (b) present the evolutions of the minimum damping ratio of the CRE for
each generation of scenario groups 1 and 2. These figures can serve as an illustration of
the convergence trajectories during the computation as well.
Figure 4.7: Convergence of the evolution
As shown in fig. 4.7, the evolutionary algorithms applied to the scenarios of group 1 and
2 converge after 230 and 426 generations respectively. The individuals which have the
best fitness are reported as follows:
• For scenario group 1,

the best individual:
800 812 848 844 828 862 800 1 1 0.7 1 0.8 1 1
the damping ratio of the CRE: 0.0785
the corresponding CRE: -0.331 ± j4.209

the best individual:
838 840 838 840 836 838 846 1 1 1 1 1 1 1
From the results it can be expected that the worst scenario appears when:
62 Chapter 4
1. the WGS are connected to the nodes, which have a relatively weak connection to
the network (normally these nodes will have a relatively big value of short circuit
impedance);
2. the WGS have high generation levels.
These conjectures serve as a direction for the study in the following step.
Step 2, choosing nodes 822, 838, 840, 848, 856, 864 and 890 for a thorough study to
seek the worst case. The reason of choosing these nodes is their relatively weak connec-
tions to the network. Further, the possible values of the wind speed are confined to 0.9
and 1.0, to give a high generation level.
Figure 4.8: Determination of the worst scenario
Fig. 4.8(a) gives the relation of the sum of the active power of 7 WGS versus the damping
ratios of the CRE for scenario group 1. In fig. 4.8(b) the relation between the generation
Chapter 4 63
level of WGS and the real parts of the CRE is shown. Fig.4.8(c) and (d) report the corre-
sponding results for scenario group 2.
Detailed information about the worst scenario is given as follows, in which, to make
the expressions concise, the genetic structure used in the previous step is adopted:

the conditions for the worst scenario:
822 822 822 822 822 822 822 1 1 1 1 1 1 1

890 890 890 890 890 890 890 1 1 1 1 1 1 1
By investigating the results, the conclusions made in the previous step for the occurrence
of the worst scenario can be improved as follows
1. all the WGS are connected to a certain node which is a weak connection point of the
network; however, this specific node varies when the operation status of the system,
such as the load level, alters;
2. the WGS are with full power generation;
3. generally when the network is lightly loaded, the situation is worse.
These results are verified in the following step.
Step 3, computing a large amount of randomly generated scenarios to verify the results
given above. In this step, for each scenario group of 1 and 2, 50000 scenarios are ran-
domly generated. Their minimum damping ratios and the corresponding generation levels
are reported in fig. 4.9(a) and (b) respectively. For convenience of the comparison, the
worst case obtained in step 1 and 2 are also marked in it respectively with the symbols of
2 and 4.
64 Chapter 4
Figure 4.9: Result verification by a large amount of random scenarios
By considering the results obtained in step 1, 2 and 3, simultaneously observing the logi-
cal relationship among them, it is reasonable to conclude that step 2 figures out the worst
scenario among all the potential candidates.
In [22], it is pointed out that for power systems, the minimum acceptable damping ratio
should be not less than 0.03. This conclusion is summarized mainly based on the experi-
ence of transmission system operations. The application of it to distribution systems may
need further verification. However, before more experiences and proper conclusions are
available, it is meaningful to take it as a reference. With the help of it, a statement can be
generalized as:
Normally the high penetration level of DG will not lead to small signal instability prob-
lems in distribution systems. Compared to the small signal stability, the steady state
voltage rise of the system is a more important factor to determine the penetration level of
DG.
Simulation Results for the Scenarios of Group 3 and 4
Studies on the scenarios of group 3 and 4 given in table 4.1 have been carried on with a
same routine as that for group 1 and 2. Summaries of the results are given as follows:
• The maximum number of the WGS allowed to connect to the network is 6.

Chapter 4 65

822 822 822 822 822 822 1 1 1 1 1 1

890 890 890 890 890 890 1 1 1 1 1 1
Comparing these results to those of scenario group 1 and 2 leads to conclusions as:
• the weak connection between the distribution and transmission system will reduce
the allowable penetration level of DG;
• the system does not have any small signal instability problem.
4.3.2 Strategies to Increase the Penetration Level of DG and the Cor-

responding Impacts on the Small Signal Stability of the SNI
Methods to Counteract the Steady State Voltage Rise
It is well recognized that the system steady state voltage rise plays an important role on
limiting the penetration level of DG. Thus many methods are proposed to cope with this
voltage rise effect [67][68][66]:
• if available, tuning the OLTC at the CCP of the transmission and distribution sys-
tem;
• adopting shunt compensating devices to adjust the reactive power injection;
• when necessary, curtailing the power generation of DG;
• implementing load control;
• selecting proper connection point for DG;
• changing line impedance.
66 Chapter 4
Curtailment of DG power generation does not meet the aim of this work; load control
may not be feasible when suitable loads, such as thermal storage heaters, are not avail-
able; considerations on network topology and parameters, such as the selection of DG
connecting point and changing line impedance, are apparent to be predicted with positive
contributions to the small signal stability. Thus, the focus will be put on the first two
methods, i.e., using OLTC and shunt compensators for voltage control.
Further, compared to shunt compensators, the OLTC is foreseen to play an ancillary role
in this work. The reasons are: on the one side, DG has a stochastic output, thus when its
penetration level is high, it will cause relatively large voltage fluctuations; on the other
side, the response of the OLTC is relatively slow, which can be tens even up to hundreds
of seconds. These two sides lead to a potential risk that before the proper operations of
the OLTC, the system has already subjected to an unallowable overvoltage lasted for an
unallowable period. Thus the functions of the OLTC in this work are simplified as: when
the load level is low, the OLTC makes the voltage of the CCP equal to 1.0; while when
the load level is high, the OLTC makes the voltage of the CCP equal to 1.05.
As for shunt compensating devices, STATCOM stands out from all the candidates for
its better operation characteristics, such as improved V-I relation, faster response time
and so on [69][70]. Hence it is adopted here to control the steady state voltage. The
relevant parameters are reported in appendix D.
Design of the Control Strategy of STATCOM
As discussed in [71] and [7], the so-called unity power factor (UPF) strategy does not
help much to increase the penetration level of WGS. Under some specific conditions it
may even adversely affect system operation characteristics. Thus here the STATCOM is
used to control nodal voltage directly rather than aim to correct the power factor and hence
affect the nodal voltage indirectly. However, calculations reveal that it is impracticable if
the set point of the STATCOM is fixed, say, Vtre f = 1.0. The reason is that the STATCOM
does not have enough capacity to maintain the nodal voltage fixed under all scenarios.
One apparent solution is to adopt a large capacity STATCOM. However, this may not be
economically justifiable. Another solution, which is used here and will be discussed later
Chapter 4 67
in detail, is to employ an adaptive set point adjustment scheme (ASPAS), to adjust the
reference voltage according to the operation state of the system.
In this work, it is observed that the fluctuation of the nodal voltage is mainly determined
by the power injection from WGS as well as the load level. Further, as for the power gen-
eration from the WGS, since the induction generator is employed, the injected reactive
power is highly dependent on the injected active power, thus the injected active power
can be used as an indicator for the generation level of the WGS. On the other hand, as
for the load level, it can be detected by summing up the power injected by WGS and that
from the CCP. More simulations further reveal that, if the load level keeps unchanged
while the output of WGS changes, though the sum of the reactive power from the WGS
and the CCP has a relatively large fluctuation, that of the active power keeps relatively
stable. Thus it can be used as an indicator for the load level. Therefore, it is proposed that
the set point of the STATCOM can be determined according to both the active power fed
by the WGS and the sum of the active power from the WGS and the CCP.
The implementation of the ASPAS works as follows: assuming a rough relation is avail-
able, which maps the active power of WGS and the sum of the corresponding active power
to an expected STATCOM set point, then if the operation state of the system changes, a
rough value of the reference voltage can be determined by referring to the relation, which
later will be called as a map. Based on it, if the operation state of the STATCOM does not
meet the prescribed criterion, the set point will further be stepwise tuned.
Map of the Active Power by the WGS and the Sum of the Active Power from the
WGS and the CCP to the set point of the STATCOM
The map consists of 2 components:
• A basic frame of the map —- The frame is actually a set of the precalculated solu-
tions of the STATCOM set point for some typical scenarios.
• Rules to get the solution of the set point for real cases —- In this work, when the
real operation state of the system is available, linear interpolations are used to get
the solution from the above mentioned basic frame.
68 Chapter 4
To form the basic frame of the map, the following factors/assumptions are considered:
• In the precalculated scenarios, two parameters vary from case to case. They are: 1.
the load level, which has options as being a minimum load demand or a maximum
load demand; 2. the wind speed, which can have a value from the set [0 0.4 0.5 0.6
0.7 0.8 0.9 1.0]. Thus there are totally 16 precalculated scenarios.
• The value of the STATCOM set point is assumed to have a resolution as 0.01 p.u.,
i.e., the value of the set point can be, for example, 0.95 or 0.96, but can not be 0.955.
• Under steady state, the STATCOM is allowed at most to provide/absorb a quantity

of reactive power equal to 70% of its capacity. This guarantees that generally the
STATCOM will have at least 30% of the reactive power capacity as a reservation,
which may be very helpful for dealing with some emergent situations.
• The STATCOM set point of the precalculated scenarios is determined by the policy
of maximizing the reservative capacity of the STATCOM.
After the mapping, the obtained solution of the set point will be further stepwise tuned if
necessary. The assumptions during the tuning are:
• The discrete tuning step is 0.01 p.u., however, the set point can not be lower than
0.9 p.u. or higher than 1.1 p.u..
• The STATCOM must have at least 30% reservative capacity under the steady state.
Results and Analyses
In order to make the study concise, by referring to previous results, it is assumed that all
the WGS and the STATCOM are connected to node 890. Simulations show that for the
scenarios of group 1 and 2 given in table 4.1, totally 9 WGS can be connected; while for
the scenarios of group 3 and 4, the corresponding penetration level is 8 WGS.
The above mentioned map for determining the set point of the STATCOM is given in
fig. 4.10(a) and (b) for group 1,2 and 3,4 respectively.
More details about fig. 4.10 are given as:

Chapter 4 69
Figure 4.10: Strategy for STATCOM set point
• The basic frame are a set of the solutions of the set point marked by markers ’+’
and ’x’, where, markers ’+’ designate the STATCOM set point when the load level
is low while markers ’x’ specify the set point when the load level is high.
• The active power shown in the figure and later discussion is scaled by the total
capacity of the installed WGS, which in this part for scenario group 1 and 2 is 18
MVA while for scenario group 3 and 4 is 16 MVA.
• The implementation of the rules to get the solution of the set point for a real case
is explained via an example given in fig. 4.10(a). Assuming an actual system
operation state is known as: 1. the active power from the WGS is 0.5; 2. the sum
of the active power from the WGS and the CCP is 0.75. As previously discussed,
the sum of the active power from the WGS and the CCP is the indicator of the load
level. Thus condition 2 can be rewritten as the load level is about 0.75. According
to the first condition, points Vs1 and Vs2 are obtained by linear interpolation. For
Vs1, it corresponds to the operation state as the active power from the WGS is 0.5,
while the load level is 1.0 (a maximum load demand). For Vs2, it corresponds to
the operation state as the active power from the WGS is 0.5, while the load level
is 0.2 (a minimum load demand). By linear interpolation according to the second
condition, the expected STATCOM set point Vs is obtained from Vs1 and Vs2,
which has a value about 1.03.
70 Chapter 4
Other relevant simulation results are:
• for scenario group 1,

the conditions for the worst scenario: the WGS produces full power generation; the
voltage of the STATCOM set point equals 1.05;
the damping ratio of the CRE: 0.0757;
the corresponding CRE: -0.312 ± j4.116.



The results demonstrate that the installation of the STATCOM increases the penetration
level of WGS. Meanwhile, as for the aspect of the small signal stability, it has been im-
proved when using the CRE and their damping ratio as an evaluation criterion.
4.3.3 Relevant Results for the SNII
When investigating the impacts of WGS on the small signal stability of the SNII, 3 factors
are expected to have a relatively large influence on the results, i.e., load level, WGS pen-
etration level and WGS generation level. The considerations for them are similar as those
Chapter 4 71
in the previous part. The main difference is about the load level, for which, a low load
level here corresponds to the load values in appendix C, while a high load level means the
magnitude of the load is 5 times bigger than that of the low load level.
Similar as in the previous part, according to the criterion of the steady state nodal voltage,
the maximum number of WGS allowed to connect to the SNII is determined as 3. It is
worth noting that in the relevant analyses, both node 1 and 2 of the SNII are modeled
as slack buses, and the difference of the voltage angles between them is determined by a
preliminary power flow calculation under the assumptions as:
• without DG
• with a low load level
• the active power injection from node 1 is 1.3 MW
• the magnitude of the voltages at node 1 and 2 is 1.05 p.u..
As for the small signal stability, the worst scenarios and their corresponding conditions
are given in the following:
• When the load level is low,

13 13 13 1 1 1
• When the load level is high,

13 13 13 1 1 1
When a STATCOM is installed, which together with the WGS is assumed to be connected
to node 13, the following results are obtained:
• The maximum number of WGS allowed to connect to the system is 4

72 Chapter 4
• When the load level is low,

• When the load level is high,

From the results it can be seen that generally the integration of WGS will not result in
small signal instability problems. The installation of STATCOM helps to increase the
penetration level of WGS as well as improve the small signal stability of the whole system.
These conform to the conclusions obtained in the previous part.
4.3.4 Conclusions and Comments

This chapter focuses on studying the impacts of DG on the small signal stability of distri-
bution systems. The simulations lead to the following conclusions:
• A high penetration level of DG will not lead to small signal instability problems.
Compared to the small signal stability, the steady state voltage rise of the system is
a more important factor to determine the penetration level of DG.
• The worst case of the small signal stability happens when the following conditions
are met: 1. all the DG are connected to a node which has a weak connection to the
network; 2. the power generation level of the DG is high; 3. the load level is low.
• Suitable reactive compensations are helpful for the integration of DG. In this work,
the application of STATCOM increases the penetration level of DG while at the
same time improves the small signal stability of the whole system.
So far, some instructive conclusions are obtained. However, some limitations about the
results exist as well. They are:
Chapter 4 73
• The results obtained in this chapter are based on the assumption that the system is
three phase symmetric and balanced. These are common assumptions when adopt-
ing the modal analysis to investigate the small signal stability. In case the unsym-
metric structure or the unbalanced operation is too heavy to be omitted, some other
more sophisticated theories have to be adopted.
• In this work, the DG are confined to the constant speed wind generation systems. To
get more general conclusions, an apparent solution is developing a universal model
of DG to execute the corresponding analyses. However, the development of such a
model will be a demanding task, not only because of the various characteristics of
different types of DG, but also the lack of practical data. The deficiency of the nec-
essary data will also challenge the verification of the universal model (if available).
All in all, to develop an universal model is very helpful if not indispensable, but it
is hard to be achieved because of the limited resources at the present stage and can
constitute a core part of the future work.
• It is well known that the normal way of wind generation is to collect several wind
turbines together to form a wind farm to generate and feed electricity into the net-
work. Thus the investigations in section 4.3.1 about each WGS connected to a dif-
ferent node seems not to agree with the practical situation. Besides, assuming the
wind blown on each wind turbine has totally independent speed also does not agree
with the real situation. However, since the object of this work is to evaluate the
impacts of DG rather than only WGS on system operations in a more general sense,
to allow the WGS to be connected to the network arbitrarily and at the same time
have totally independent power generation levels (equivalent to totally independent
wind speeds) may give a more suitable approximation to the potential operation of
DG in the future thus leads to some more general conclusions.
• In section 4.3.2, the operation of the OLTC has been simplified. Efforts can be
made in the future to coordinate the OLTC and the STATCOM more effectively to
increase the penetration level of DG as well as improve the system operations.
Chapter 5
Impacts of DG on the Voltage Stability

of Distribution Systems
Voltage stability refers to the ability of a power system to maintain steady voltages at all
buses in the system after being subjected to a disturbance from a given initial operating
condition. It is reasonable to expect that the installation of DG will enhance the voltage
stability of a distribution system. In this chapter, this conjecture will be verified with the
method of continuation power flow to investigate the voltage stability margin improve-
ment after the integration of DG. Further, the effectiveness of some possible operation
modes of DG on the improvement of the voltage stability is discussed.
In this chapter, section 5.1 presents a detailed description of the implementation of the
continuation method. Section 5.2 gives information about the assumptions and models
for the analyses in this chapter. Simulations, results and comments are available in sec-
tion 5.3.
The work of this chapter shows that the integration of DG is helpful to improve the volt-
age stability of distribution systems. Further, DG operation mode has a great influence on
the improvement of the voltage stability. The greatest improvement is obtained when the
DG is operated under the constant nodal voltage mode.
Chapter 5 75
5.1 Implementation of the Continuation Method

Continuation methods are important tools for implementation of branch tracing or path
following. A most frequently and successfully used continuation method may be the one
based on a so-called ’predictor - corrector’ procedure. [72] well documents its applica-
tion on general nonlinear systems. While [73] [74] [75] thoroughly describe the way of
its implementation in power system analyses. In the following, its principles are reviewed
in detail by taking power systems as an example.
A quasi steady power system can be modeled as
0 = f (x, λ) (5.1)
where, x is a n-dimensional vector, and depending on the modeling level, x can only in-
clude nodal voltages and corresponding angles, or further consist of some state variables;
while λ is normally a scalar variable to be used to control a continuative variation of cer-
tain parameters.
For example, when investigating the power system loadability on a peculiar direction,
the system load can be modeled as:
PL = PL0 + λPd
QL = QL0 + λQd (5.2)
where PL0 and QL0 are the original load of the system; while Pd and Qd designate a par-
ticular direction of load variation.
In order to find a solution of eqn. (5.1), firstly a guess of the solution must be made
(this actually is done by the predictor), then from the guess (the initial value) the real
solution has to be found (this is fulfilled by the corrector). The following parts give the
details of this procedure.
Predictor
Predictors normally can be divided into two classes:
76 Chapter 5
• ODE-based method —- This method exploits the current solution and the concerned
derivatives to form the initial guess of the next solution. When only the first order
derivative is used, it forms the so-called tangent predictor.
• Polynomial extrapolation based method —- This method uses the current as well as
the previous solutions to predict the next solution. When a zero order polynomial is
adopted, this method equals to take the present solution directly as the initial guess
of the next solution. When a first order polynomial (a line) is adopted, it leads to
the so-called secant predictor.
Tangent Predictor
Differentiating eqn. (5.1) leads to
0 = fx dx + fλ dλ (5.3)
where fx and fλ are the corresponding Jacobian matrices.
When fx is nonsingular, eqn. (5.3) leads to
dx
= −( fx )−1 fλ (5.4)
dλ
Note that eqn. (5.4) gives a increment direction of x when λ increases. Assuming the
present solution as (x0 , λ0 ), given a stepsize as ∆λ, thus the increment of x is
∆x = −( fx |x=x0 ,λ=λ0 )−1 fλ |x=x0 ,λ=λ0 ∆λ (5.5)
and therefore the guess of the next solution can be predicted as
p
x1 = x0 + ∆x
p
λ1 = λ0 + ∆λ (5.6)
where, the superscript p means a predicted value.
However, this procedure fails when the bifurcation point is reached, where the Jacobian
matrix fx is singular. This problem can be solved by switching λ in the previous equations
to an alternative parameter. It can be achieved as follows:
Chapter 5 77
When the aforesaid procedure approaches the bifurcation point, it is noticed that when
λ has a small increment, some other parameters may subject to a sharp change. Under
this situation, the one which has the largest relative increment, say xi , can be used to sub-
stitute λ to implement the aforesaid procedure in a similar way.
Eqn. (5.3) will be rearranged as
0 = fx̂ d x̂ + fxi dxi (5.7)
where, x̂ = [x1 ... xi−1 xi+1 ... xn λ ]T .
Thus the incremental direction of x̂ with respect to xi is

d x̂
= −( fx̂ )−1 fxi (5.8)
dxi
Assuming the current solution as (x̂0 , xi0 ), given a stepsize as ∆xi , hence the corresponding
increment of x̂ is
∆x̂ = −( fx̂ |x̂=x̂0 ,xi =xi0 )−1 fxi |x̂=x̂0 ,xi =xi0 ∆xi (5.9)
and therefore the next point can be predicted as

p
x̂1 = x̂0 + ∆x̂
p
xi 1 = xi 0 + ∆xi (5.10)
When xi is properly selected, the Jacobian matrix fx̂ is guaranteed to be nonsingular even
at the bifurcation point, thus the feasibility of eqn. (5.9).
Secant Predictor
Assuming the present and previous solutions of eqn. (5.1) are known as (x1 , λ1 ) and (x0 ,
λ0 ), then the secant predictor adopts the first order polynomial to predict the next point
as:
p
x2 = x1 + h(x1 − x0 )
p
λ2 = λ1 + h(λ1 − λ0 ) (5.11)
78 Chapter 5
where, h is a proper stepsize.
Corrector
A corrector is actually an algorithm to solve the real solution (x, λ) of eqn. (5.1) with the
available initial guess of (x p , λ p ), which is from the predictor. Newton method is popu-
larly utilized in various correctors for its effectiveness of solving nonlinear equations.
Observation of eqn. (5.1) shows that it includes n equations but n + 1 variables. This
means another extra equation is necessary for the solution. The procedure of forming an-
other required equation is called parameterization. The normally used parameterization
methods include: 1. local parameterization; 2. perpendicular intersection parameteriza-
tion; 3. arclength parameterization. In the following they are explained in detail.
Local Parameterization
In local parameterization, either one element of the n−dimensional vector x or the scalar
λ is forced to be a fixed value. That is, one of the following equations will be added
together with eqn. (5.1) for seeking a practical solution.
λ − λp = 0 (5.12)
p
xi − xi = 0 (5.13)
Normally eqn. (5.12) is adopted at the beginning of the procedure. With the approaching
of the bifurcation point, eqn. (5.13) will be used instead.
Perpendicular Intersection Parameterization
It is assumed that in perpendicular intersection parameterization, for the current solution

p p
(x0 , λ0 ), the prediction of the next solution (x1 , λ1 ), and the real value of the next solution
p p
(x1 , λ1 ), the hyperplane determined by (x0 , λ0 ) and (x1 , λ1 ) is perpendicular to the one
p p
determined by (x1 , λ1 ) and (x1 , λ1 ). That is, the following condition is a complement of
eqn. (5.1).
p p p p
[x1 − x0 ]T ∗ [x1 − x1 ] + (λ1 − λ0 ) ∗ (λ1 − λ1 ) = 0 (5.14)
Chapter 5 79
Arclength Parameterization
The arclength parameterization method introduces the arclength s on the solution curve
as a new parameter. From the current solution (x0 , λ0 ) to the next solution (x1 , λ1 ), the
corresponding increment of the arclength, i.e., ∆s, is under control. This leads to the
following additional equation for eqn. (5.1).
n
∑ (xi1 − xi0)2 + (λ1 − λ0)2 − (∆s)2 = 0 (5.15)
i=1
Step Length Control

One apparent step length control strategy is to take a constant step length, say, with
∆λ = 0.1 during the procedure. If the step length is small enough, it may work suc-
cessfully for a wide range of problems. However this will often lead to inefficient com-
putation, since too many steps are spent on the ’flat’ part of the solution curve. On the
contrary, if a larger step length is adopted, though the time for branch tracking on the flat
part may be reduced, when the curvature increases, it can cause the predicted value (from
the predictor) deviating far from the true solution, thus the time for the corrector computa-
tion increases. In some extreme case, the corrector may even be subjected to divergence.
All these imply that an adaptive step length control strategy has to be used.
In [72], it is pointed out that ”the total amount of work of a continuation procedure de-
pends on the average step length in a somewhat convex manner: The continuation is
expensive both for very short steps (too many steps) and for very large steps (slow or no
convergence of the corrector). Its costs are moderate for a certain medium step length,
which is related to an optimal number Nopt of iterations of a corrector” (quoted from
[72]). Thus the relation of the predetermined Nopt and the actual iteration number N j of
the present corrector computation can be adopted to adjust the step length for the next
computation adaptively. When Nopt > N j , the step length is considered to be increased;
while when Nopt < N j , it will be reduced .
Another option of step length control strategy may take the curvature of the solution
curve into account as well. This can be done by assuming λ is the controlled continuative
parameter, thus the scaled directions of x with respect to λ at the previous solution of
80 Chapter 5
(x0 , λ0 ) and the current solution of (x1 , λ1 ) can be used to determine the step length as
hc
h= (5.16)
||J1 − J0 ||
where J is the scaled direction, hc a constant mainly depending on experience, h the step
length, which in practice may be bounded, say, as

 0.1 if h < 0.1


h= h if 0.1 ≤ h ≤ 1 (5.17)


1 if h > 1

Eqn. (5.16) reveals that at the flat part of the solution curve, since the directions will not
change so much, their difference thus will be relatively small, therefore a large step length
is obtained. However, when the curve turns sharply, the distance between two scaled
directions increases quickly as well, thus a small step length is accordingly adopted.
5.2 Models of Power System Components

In this chapter, the steady state models are used:
• The transmission system is modeled as an infinite bus.
• The electrical network is mathematically modeled as a nodal admittance matrix.
• The load is modeled as a constant PQ component.
• The DG is modeled as a PQ generator. However, when enough capacities of the

shunt compensator are assumed in operation, they together will be modeled as a PV
generator.
The DG is modeled as a PQ generator because normally it is not controlled by system

operators. This model may be suitable for those types of DG which have relatively stable
performance, such as fuel cell generations and combined heat and power plants. But for
other DG, such as WGS and photovoltaic generation systems, it is not true. However, with
the development of various feasible energy storage systems, the output of the aforesaid
DG may be expected to be well stabilized in the future, thus makes the PQ model to be
acceptable.
Chapter 5 81
Impacts of DG on the Voltage Stability of the SNI

Based on the discussion of the previous chapter, the following assumptions are adopted
to make the analyses concise as well as meaningful.
• The voltage at the node 800 is fixed as 1.05 p.u..
• The DG is connected to node 890, which has an output of the active power as
12MW. However, the corresponding reactive power is determined according to
three options of the power factor: unitary power factor, capacitive power factor
0.95 and inductive power factor 0.95.
• The parameters of the load given in appendix B are adopted here as the load varia-
tion direction. Hence in the following analysis, the load is changed as:
PL = λPL0
QL = λQL0 (5.18)
where λ is a positive real number. PL0 and QL0 are the initial conditions of the load
given in the precious chapter.
After the continuation power flow calculation, the relation between the nodal voltage of
bus 890 and the total active power of the load is presented in fig. 5.1, which reveals
that generally the integration of DG is helpful for the voltage stability of the distribution
system. Further, the largest gain of the voltage stability margin is obtained when the con-
nected DG is operated under the capacitive power factor mode.
If enough shunt compensators are installed for nodal reactive power injection control,
then they may offer another optional operation mode – keep the concerned nodal volt-
age constant. Compared to the capacitive power factor mode, this mode can make the
voltage stability properties even better. Fig. 5.2 shows the results using the relation of
the nodal voltage of bus 888 and the total active power of the load as an example, where
the nodal voltage 890 is assumed to be as 1.05 p.u.. Besides, the corresponding curve of
the aforementioned capacitive power factor mode is also drawn in it for convenience of
comparison.
82 Chapter 5
Figure 5.1: Loadability of the SNI for DG under different operation modes
Figure 5.2: Comparison of the loadability of the SNI for DG under the operation modes
of constant voltage and capacitive power factor
Chapter 5 83
Impacts of DG on the Voltage Stability of the SNII
Similar as in the previous part, some assumptions are adopted here to simplify the analy-
ses:
• The voltage at the node 1 and 2 are fixed as 1.05 p.u., and the difference of their
arguments is also fixed. Further, when enough shunt compensation is available, the
concerned nodal voltage is assumed to be 1.05 p.u..
• The DG as well as the shunt compensator (if available) is connected to node 13.
The active power output of the DG is 6MW. There are 4 options of the operation
mode: unitary power factor, capacitive power factor 0.95, inductive power factor
0.95 and constant voltage.
• The parameters of the load given in appendix C are the load variation direction.
Fig. 5.3 gives the results of the continuation power flow analyses for the scenarios of
without DG and DG as a PQ generator. Since the SNII is fed of electricity by the trans-
mission system simultaneously at two nodes, its voltage stability is much better compared
to that of the SNI. In fig. 5.3, only part of the whole P-V curve for node 13 is given. To
compute a complete P-V curve is not necessary, since the bifurcation point of the SNII
will appear when the load is so heavy as to override the thermal limit of the cable. Thus to
determine such a bifurcation point is physically meaningless. Without information about
the location of the bifurcation point, an obvious observation about the voltage stability
improvement is not available. However, fig. 5.3 shows that the nodal voltage of bus 13 is
increased after the integration of the DG, hence it is reasonable to expect that the voltage
stability of the whole system has been improved.
As for the scenarios of the DG modeled as a PV generator, the results are given in fig.
5.4. It indicates that compared to the mode of capacitive power factor, when the DG is
operated under the constant voltage mode, the nodal voltage drop will be smaller as the
increment of the load level. Thus it is reasonable to predict that the constant voltage mode
will be more helpful for the voltage stability of the system.
84 Chapter 5
Figure 5.3: Loadability of the SNII for DG under different operation modes
Figure 5.4: Comparison of the loadability of the SNII for DG under the operation modes
of constant voltage and capacitive power factor
Chapter 5 85
Comments
So far, a positive expectation of DG to improve the voltage stability is given. However,
this work may be improved when the following aspects are taken into account:
• As in the previous chapter, the results obtained in this chapter are also based on a
three phase symmetric system under balanced operation. In case of investigations
under unsymmetric structure or unbalanced operation, methods proposed in [76]
are applicable.
• Compared to traditional electricity generation, the operation of DG contains more

stochastic characteristics. Thus to investigate their impacts on the voltage stability
from a more stochastic standpoint may be more suitable and helpful. However,
research in this aspect is still open since the lack of practical experience and data.
Chapter 6
Conclusions and Future Works
6.1 Conclusions
Concerning the studies of the eigenvalues of the state space model of distribution
networks, the following conclusions are drawn:
• There is a close relation between the eigenvalues of the network under abc and dq0
reference frame. Assuming a symmetric network has 3N state variables, thus it
has 3N eigenvalues. Among them, N pairs of the complex-conjugate/double-real
eigenvalues under abc frame, relate only to the network dq components under dq0
frame. These eigenvalues under dq0 frame can be gotten by shifting the counter-
parts under abc frame with frequency of ±ω0 . The remaining N eigenvalues will
keep unchanged under both of the two frames.
• The parameters of the lines/cables of the electrical network are the predominant
factors to determine the eigenvalues of the network. Although the network topology
has an influence on the eigenvalues, the parameters of the lines/cables are of much
greater influence.
• The distribution of the eigenvalues of the network is bounded by the relevant values
of its lines/cables. Assuming a distribution network consists of M different types
of lines, then the maximum and minimum values of the real parts of their CE span
network.
Chapter 6 87
• In a real distribution network in Europe, the eigenvalues will have the real parts
within the range from -69.3 to -33055. This is based on a thorough investigation of
the data of the lines/cables given in Neplan [46].
• The impacts of the apparently successive stochastic disturbances to the network

can be studied independently and the nodal admittance matrix is acceptable as the
model of a distribution network. The investigation on the eigenvalues of a distribu-
tion network reveals that the electrical dynamics die off very fast. Hence a steady
state assumption for modeling the distribution network is acceptable. The appar-
ently successive stochastic disturbances have nearly independent influences on the
network.
Concerning the impacts of DG on the small signal stability of distribution systems,

the following conclusions are drawn:
• A high penetration level of DG will not lead to small signal instability problems.
Compared to the small signal stability, the steady state voltage rise of the system is
a more important factor to determine the penetration level of DG.
• The worst case of the small signal stability happens when the following conditions
are met: 1. all the DG are connected to a node which has a weak connection to the
network; 2. the power generation level of the DG is high; 3. the load level is low.
• Suitable reactive compensations are helpful for the integration of DG. In this work,
the application of STATCOM increases the penetration level of DG while at the
same time improves the small signal stability of the whole system.
Concerning the impacts of DG on the voltage stability of distribution systems, the

following conclusions are drawn:
• In general the integration of DG improves the voltage stability of the distribution

system.
• The operation mode of DG has a great influence on the improvement of the volt-
age stability. The largest gain of the voltage stability margin is obtained when the
connected DG is operated under the constant nodal voltage mode.
88 Chapter 6
6.2 Future Works

Works carried out to present stage result in some useful conclusions. However, improve-
ments are expectable in various aspects. Following comments give some indications for
the future works in general.
• Development of a universal model of DG - A universal model of DG is very help-

ful if not indispensable for concerned analyses. Though presently hard to achieve,
with the accumulation of practical data and experience, it may be feasible in the
future.
• More proper considerations of the stochastic characteristics of DG - Compared

to conventional generation, DG is more stochastic in its operation. Efforts to include
this aspect more properly in concerned analyses are meaningful.
• Researches under three phase unsymmetric/unbalanced conditions - Distribu-

tion systems are apt to be subjected to unsymmetric/unbalanced conditions when
compared to transmission systems. Thus researches carried on under these condi-
tions will be more helpful for utilities. However, this is by no means an easy work.
New and appropriate methods for this aim may need to be developed.
Appendix
A. Explanation of Equation (3.13)
Observation
For networks with a radial structure, if the shunt admittance of their lines/cables is negli-
gible, then through a systematic routine their state matrix can be formed in a generalized
form as:
İ2(n−1)×1 = A2(n−1)×2(n−1) I2(n−1)×1 + B2(n−1)×2 v2×1 + C2(n−1)×2(n−1) V2(n−1)×1 (1)
where
n is the total number of the nodes of the network;
I is the set of the state variables consisting of the dq components of the deliberately se-
lected branch currents;
both v and V are sets consisting of the dq components of the specially selected nodal
voltages.
Example
The network given in fig. 1 is adopted as an example to explain the observation, in which,
node 1 is a SN, node 4 and 5 are TN, and node 2, 3, 6, 7 and 8 are EN.
For the structure enclosed by the dotted curve in fig. 1, the following equations are valid:

 L6 di6d − L i + R i
ω0 dt 6 6q 6 6d = v5d − v6d
L di (2)
6q
 6
ω0 dt + L6 i6d + R6 i6q = v5q − v6q

L7 di7d

ω0 dt − L7 i7q + R7 i7d = v5d − v7d
L7 di7q (3)

ω0 dt + L7 i7d + R7 i7q = v5q − v7q

L8 di8d

ω0 dt − L8 i8q + R8 i8d = v5d − v8d
L8 di8q (4)

ω0 dt + L8 i8d + R8 i8q = v5q − v8q

L5 di5d

ω0 dt − L5 i5q + R5 i5d = v4d − v5d
L5 di5q (5)

ω0 dt + L5 i5d + R5 i5q = v4q − v5q
90 Appendix
Figure 1: Example for explanation of the general form of state matrix of radial
impedance networks
with algebraic equations
i5d = i6d + i7d + i8d

i5q = i6q + i7q + i8q (6)
where
im (m = 5,6,7,8) are branch currents flowing from the SN to the TN, or from the TN to the
EN;
Rm and Lm are the corresponding branch impedance.
For example, i5 is the current on the branch with a direction from node 4 to 5, while
R5 and L5 the corresponding impedance.
It can be noticed that equations (2),(3) and (4) have the following general form
İi = Ai Ii + Bi Vstart + Ci Vend (7)
where     
− RLi ωi 0 ω0 ω0
Li
ω0
Li
Ai =   , Bi =   , Ci = −  
−ω0 − RLi ωi 0 ω0
Li
ω0
Li
Appendix 91
i = 6, 7, 8
and the contents of Ii , Vstart , Vend can be clarified via taking eqn. (2) as an example, in
which, Ii = [i6d i6q ]T , Vstart = [v5d v5q ]T and Vend = [v6d v6q ]T .
Then equations (2),(3) and (4) can be organized as
İ678 = A678 I678 + B678 V5 + C678 V678 (8)
where      
A6 B6 C6
     
A678 = 
 A7  , B678 =  B7  , C678 = 
    C7 

A8 B8 C8
I678 = [i6d i6q i7d i7q i8d i8q ]T
V5 = [v5d v5q ]T
V678 = [v6d v6q v7d v7q v8d v8q ]T
From eqn. (5) and (6) the following equation can be deduced
V5 = DV4 + Ei İ678 + Fi I678 (9)
where   
1 1 1 1
D=  , Ei = − L i  
ω0
1 1 1 1
   
1 1 1 1 1 1
Fi = −Ri   + Li  
1 1 1 −1 −1 −1
V4 = [v4d v4q ]T
with i = 4.
Eqn. (8) and (9) result
İ678 = Ā678 I678 + B̄678 V4 + C̄678 V678 (10)
where
Ā678 = (16×6 − B678 E4 )−1 (A678 + B678 F4 ),
92 Appendix
B̄678 = (16×6 − B678 E4 )−1 B678 D,

C̄678 = (16×6 − B678 E4 )−1 C678
with
16×6 an unit matrix with a dimension of 6 × 6, superscript -1 the inversion of the corre-
sponding matrix.
The mathematical model of the whole network shown in fig. 1 consists of eqn. (10)
as well as the following equations

 L2 di2d − L i + R i
ω0 dt 2 2q 2 2d = v4d − v2d
L di (11)
2q
ω0 dt + L2 i2d + R2 i2q = v4q − v2q
 2

L3 di3d

ω0 dt − L3 i3q + R3 i3d = v4d − v3d
L3 di3q (12)

ω0 dt + L3 i3d + R3 i3q = v4q − v3q

L4 di4d

ω0 dt − L4 i4q + R4 i4d = v1d − v4d
L4 di4q (13)

ω0 dt + L4 i4d + R4 i4q = v1q − v4q
and
i4d = i2d + i3d + i6d + i7d + i8d

i4q = i2q + i3q + i6q + i7q + i8q (14)
Equations (10),(11) and (12) can be arranged as
İ23678 = A23678 I23678 + B23678 V4 + C23678 V23678 (15)
where      
A2 B2 C2
     
A23678 = 
 A3  , B23678 =  B3  , C23678 = 
    C3 

Ā678 B̄678 C̄678
I23678 = [i2d i2q i3d i3q i6d i6q i7d i7q i8d i8q ]T
V4 = [v4d v4q ]T
V23678 = [v2d v2q v3d v3q v6d v6q v7d v7q v8d v8q ]T
Appendix 93
Combining eqn. (13) and (14) leads to
V4 = DV1 + Ei İ23678 + Fi I23678 (16)
where
V1 = [v1d v1q ]T
D, Ei , Fi the same as in eqn. (9) with i = 1.
Eqn. (15) and (16) result in
İ23678 = Ā23678 I23678 + B̄23678 V1 + C̄23678 V23678 (17)
where
Ā23678 = (110×10 − B23678 E1 )−1 (A23678 + B23678 F1 ),
B̄23678 = (110×10 − B23678 E1 )−1 B23678 D,
C̄23678 = (110×10 − B23678 E1 )−1 C23678
with
110×10 an unit matrix with a dimension of 10 × 10.
94 Appendix
B. Data of the SNI
Node Length(km) Node Length(km) Node Length(km)

800 - 802 0.7864 802 - 806 0.5273 806 - 808 9.8238
808 - 810 1.7691 808 - 812 11.4301 812 - 814 9.0618
814 - 850 0.0030 816 - 818 0.5212 816 - 824 3.1120
818 - 820 14.6763 820 - 822 4.1880 824 - 826 0.9236
824 - 828 0.2560 828 - 830 6.2302 830 - 854 0.1585
832 - 858 1.4935 834 - 860 0.6157 834 - 842 0.0853
836 - 840 0.2621 836 - 862 0.0853 842 - 844 0.4115
844 - 846 1.1095 846 - 848 0.1615 850 - 816 0.0945
852 - 832 0.0030 854 - 856 7.1111 854 - 852 11.225
858 - 864 0.4938 858 - 834 1.7770 860 - 836 0.8169
862 - 838 1.4813 888 - 890 3.2187 832 - 888 1.6094
Table 1: Line Segment Data
R = 0.128 Ohm/km, X = 0.122 Ohm/km, B = 116.239 uS/km

R0 = 1 Ohm/km, X0 = 1 Ohm/km, B0 = 116.239 uS/km
Table 2: Parameters of the Line Conductor†
Node Ph-A Ph-B Ph-C Node Ph-A Ph-B Ph-C

kVar kVar kVar kVar kVar kVar
844 100 100 100 848 150 150 150
Table 3: Shunt Capacitors
† The parameters are got from a cable made in Germany, the type of which is ’20kV NEKEBA 3×150’.
Appendix 95
Node P Q Node P Q Node P Q Node P Q

kW kVar kW kVar kW kVar kW kVar
802 275 145 806 275 145 808 80 40 810 80 40
816 25 10 818 170 85 820 845 435 822 675 350
824 245 120 826 200 100 828 55 25 830 485 215
832 75 35 834 890 450 836 610 215 838 140 70
840 470 310 842 45 25 844 4320 3290 846 340 170
848 715 535 854 20 10 856 20 10 858 245 125
860 1740 1060 862 140 70 864 10 5 890 2250 1125
Table 4: Loads
96 Appendix
C. Data of the SNII
Node Type of the cable Length(km) Node Type of the cable Length(km)
1-3 NAKBA 3×150 0.01 19 - 20 NA2XS(F)2Y 3×1×185 0.9086
3-4 NAKBA 3×150 0.2409 20 - 21 NAKBA 3×240 0.2885
4-5 NAKBA 3×150 0.226 21 - 22 NAKBA 3×240 0.8286
1-6 NAKBA 3×150 0.217 22 - 23 NA2XS(F)2Y 3×1×185 0.3374
6-7 NKBA 3×70 0.142 15 - 24 NA2XS(F)2Y 3×1×185 0.7855
7-8 NKBA 3×70 0.1334 24 - 25 NA2XS(F)2Y 3×1×185 0.9508
8-9 NA2XS(F)2Y 3×1×185 0.0416 25 - 26 NAKBA 3×150 0.4218
9 - 10 NKBA 3×70 0.2549 26 - 27 NAKBA 3×150 0.0982
10 - 11 NKBA 3×70 0.4249 2 - 28 NA2XS(F)2Y 3×1×185 0.4994
11 - 12 NAKBA 3×70 0.3664 28 - 19 NA2XS(F)2Y 3×1×185 0.9216
12 - 13 NAKBA 3×70 0.6168 19 - 29 NA2XS(F)2Y 3×1×185 0.4956
1 - 14 NA2XS(F)2Y 3×1×185 0.2619 29 - 30 NA2XS(F)2Y 3×1×185 0.2
14 - 15 NA2XS(F)2Y 3×1×185 0.4069 30 - 31 NA2XS(F)2Y 3×1×185 0.2
15 - 16 NAKBA 3×150 0.685 20 - 32 NAKBA 3×150 0.0661
16 - 17 NAKBA 3×150 1.6537 20 - 33 NAKBA 3×150 0.2606
17 - 18 NAKBA 3×150 0.71 33 - 34 NAKBA 3×150 0.2034
18 - 19 NAKBA 3×95 0.3076
Table 1: Line Segment Data

Appendix 97
Type of the cable R (Ohm/km) X (Ohm/km) B (uS/km)

NAKBA 3×70 0.446 0.104 103.673
NAKBA 3×95 0.323 0.1 116.239
NAKBA 3×150 0.21 0.094 141.372
NAKBA 3×240 0.13 0.089 169.646
NA2XS(F)2Y 3×1×185 0.164 0.161 127.235
NKBA 3×70 0.271 0.104 103.673
Table 2: Parameters of the Cable
Node S (kVA) Node S (kVA) Node S (kVA) Node S (kVA) Node S (kVA)
3 150 9 30 14 210 25 40 30 30
4 80 10 40 21 30 26 90 31 50
5 120 11 200 22 80 27 50 32 50
7 70 12 20 23 20 28 100 33 100
8 50 13 60 24 60 29 100 34 70
Table 3: Loads†
† All the loads are assumed to have a power factor of 0.95.

98 Appendix
D. Data for System Components
Variable Description Value Variable Description Value

- rated power 400 MVA - rated voltage 25 KV
ra armature resistance 0.001 p.u. xl leakage reactance 0.05 p.u.
xd d-axis synchronous reactance 1.9 p.u. M inertia constant 10 s
xd0 d-axis transient reactance 0.302 p.u. D damping coefficient 10 p.u.
Table 1: Parameters of the Equivalent Synchronous Generator

- rated power 2 MVA - rated voltage 25 KV
- rated wind speed 15 m/s ρ air density 1.225 kg/m3
rS stator resistance 0.01 p.u. xR stator reactance 0.1 p.u.
rR rotor resistance 0.01 p.u. xR rotor reactance 0.08 p.u.
xm magnetizing reactance 3 p.u. Hwr wind turbine inertia 2.5 s
Hm rotor inertia 0.5 s Ks shaft stiffness 0.3 p.u.
p number of poles 4 ηGB gear box ratio 1:89
R blade length 37.5 m bc compensating capacitor 0.5 p.u.
Table 2: Parameters of the Constant Speed Wind Turbine Generation System

Appendix 99

- rated power 5 MVar - rated AC voltage 2 KV
- rated DC voltage 4 KV Rt AC circuit resistance 0.01 p.u.
Xt AC circuit reactance 0.1 p.u. Rdc DC circuit resistance 400 p.u.
Cdc DC circuit capacitance 0.25133 p.u. Kacp - 10
Kaci - 10 Kdcp - 10
Kdci - 10
Table 3: Parameters of STATCOM

Bibliography
[1] R.A.Schlueter, G.L.Park, M.Lotfalian, et al, ’Modification of Power System Oper-

ation for Significant Wind Generation Penetration’, IEEE Transactions on Power
Apparatus and Systems, Vol. PAS-102, No.1, pp. 153-161, January 1983.
[2] A. Bhowmik, A. Maitra, S.M. Halpin, J.E. Schatz, ’Determination of allowable pen-
etration levels of distributed generation resources based on harmonic limit consider-
ations’, IEEE Transactions on Power Delivery, Vol. 18, pp. 619-624, April 2003.
[3] W. Freitas, J.C.M. Vieira, A. Morelato, W. Xu, ’Influence of excitation system con-
trol modes on the allowable penetration level of distributed synchronous generators’,
IEEE Transactions on Energy Conversion, Vol. 20, pp. 474-480, June 2005.
[4] J.G. Slootweg, ’Wind Power Modelling and Impact on Power System Dynamics’,
PhD thesis, Technische Universiteit Delft, Netherlands, 2003
[5] R.T. Guttromson, ’Modeling distributed energy resource dynamics on the trans-
mission system’, IEEE Transactions on Power Systems, Vol. 17, pp. 1148-1153,
November 2002.
[6] M.K. Donnelly, J.E. Dagle, D.J. Trudnowski, G.J. Rogers, ’Impact of the distributed
utility on transmission system stability’, IEEE Transactions on Power System, Vol.
11, pp. 741-746, May 1996.
[7] W. Freitas, A. Morelato, Wilsun Xu, F. Sato, ’Impacts of AC Generators and DSTAT-
COM devices on the dynamic performance of distribution systems’ , IEEE Transac-
tions on Power Delivery, Vol. 20, Part 2, pp. 1493-1501, April 2005.
[8] W. Freitas, J.C.M. Vieira, A. Morelato, L.C.P. daSilva, V.F. da Costa, F.A.B. Lemos,
’Comparative Analysis Between Synchronous and Induction Machines for Dis-
Bibliography 101
tributed Generation Applications’, IEEE Transactions on Power Systems, Vol. 21,

pp. 301-311, February 2006.
[9] W. Freitas, L.C.P. Da Silva, A. Morelato, ’Small-Disturbance Voltage Stability of

Distribution Systems With Induction Generators’, IEEE Transactions on Power Sys-
tems, Vol. 20, pp. 1653-1654, August 2005.
[10] S. Jang, K. Kim, ’An islanding detection method for distributed generations using
voltage unbalance and total harmonic distortion of current’, IEEE Transactions on
Power Delivery, Vol. 19, pp. 745-752, April 2004.
[11] Z. Ye, A. Kolwalkar, Y. Zhang, P. Du, R. Walling, ’Evaluation of anti-islanding

schemes based on nondetection zone concept’, IEEE Transactions on Power Elec-
tronics, Vol. 19, pp. 1171-1176, September 2004.
[12] F. Katiraei, M.R. Iravani, P.W. Lehn, ’Micro-grid autonomous operation during and
subsequent to islanding process’, IEEE Transactions on Power Delivery, Vol. 20, pp.
248-257, January 2005.
[13] N. Hatziargyriou, E. Karapidakis, and D. Hatzifotis, ’Frequency stability of power

system in large islands with high wind power penetration’, Bulk Power Syst. Dy-
namics Control Symp.-IV Restructuring, Vol. PAS-102, Santorini, Greece, pp. 24-
28, August 1998.
[14] Y. Mao, K.N. Miu, ’Switch placement to improve system reliability for radial distri-
bution systems with distributed generation’, IEEE Transactions on Power Systems,
Vol. 18, pp. 1346-1352, November 2003.
[15] N. Jenkins, R. Allan, P. Crossley, D. Kirschen and G. Strbac, ’Embedded Genera-

tion’, The Institution of Electrical Engineers, London, United Kingdom, 2000.
[16] N. Hadjsaid, J.F. Canard, F. Dumas, ’Dispersed generation impact on distribution

networks’, IEEE Transactions on Computer Applications in Power, Vol. 12, pp. 22-
28, April 1999.
[17] ’Assessment of Distributed Generation Technology Applications’, available at

http://www.distributed-generation.com/library.htm, visited on January 2006.
102 Bibliography
[18] R.C. Dugan, S.K. Price, ’Issues for distributed generation in the US’ IEEE Power
Engineering Society Winter Meeting, 2002.
[19] J. Machowski, J.W. Bialek, J.R. Bumby, ’Power System Dynamics and Stability’,
John Wiley & Sons Ltd, 1997.
[20] P. Kundur, J. Paserba, V. Ajjarapu, G. Andersson, A. Bose, C. Canizares, N. Hatziar-

gyriou, D. Hill, A. Stankovic, C. Taylor, T. Van Cutsem, V. Vittal, ’Definition and
classification of power system stability, IEEE/CIGRE joint task force on stability
terms and definitions’, IEEE Transactions on Power Systems, Vol. 19, pp. 1387-
1401, August 2004.
[21] P. Kundur, ’Power System Stability and Control’, McGraw-Hill, New York, 1994.
[22] ’Analysis and Control of Power System Oscillations’, CIGRE Report of Task Force
07, Advisory Group 01, Study Committee 38, December 1996.
[23] E. Handschin, Skriptum zur Vorlesung ’Elektrische Energietechnik 1/2’, Dortmund

University, WS 1998/1999.
[24] P. Kundur, L. Wang, ’Small signal stability analysis: experiences, achievements, and
challenges’, IEEE International Conference on Power System Technology, 2002.
[25] Y.V. Makarov, Z.Y. Dong, D.J. Hill, ’A general method for small signal stability
analysis’, IEEE Transactions on Power Systems, Vol. 13, pp. 979-985, August 1998.
[26] S. Gomes, N. Martins, C. Portela, ’Computing small-signal stability boundaries for

large-scale power systems’, IEEE Transactions on Power Systems, Vol. 18, pp. 747-
752, May 2003.
[27] M. Pavella, P.G. Murthy, ’Transient Stability of Power Systems, Theory and Prac-
tice’, John Wiley&Sons, 1994.
[28] A.A. Fouad and V. Vittal, ’Power System Transient Stability Analysis: Using the
Transient Energy Function Method’, Prentice-Hall, 1991.
[29] T.V. Cutsem, C. Vournas, ’Voltage Stability of Electric Power Systems’, Kluwer
Academic Publishers, 1998.
Bibliography 103
[30] A. Berizzi, P. Marannino, M. Merlo, M. Pozzi, F. Zanellini, ’Steady-state and dy-

namic approaches for the evaluation of loadability margins in the presence of sec-
ondary voltage regulation’, IEEE Transactions on Power Systems, Vol. 19, pp. 1048-
1057, May 2004.
[31] Federico Milano, PSAT, http://www.power.uwaterloo.ca/∼fmilano/.
[32] F. Milano, ’An Open Source Power System Analysis Toolbox’, IEEE Transactions
on Power Systems, Vol. 20, pp. 1199-1206, August 2005.
[33] E.G. Cate , K Hemmaplardh, J.W. Manke, and D.P. Gelopulos, ’Time frame no-
tion and time response of the models in transient, mid-term and long-term stability
programs’, IEEE Transactions on Power Apparatus and Systems 103(1): 143-151,
1984.
[34] P.W. Sauer, D.J. LaGesse, S. Ahmed-Zaid, and M.A. Pai, ’Reduced Order Modeling
of Interconnected Multimachine Power Systems Using Time-Scale Decomposition’,
IEEE Transactions on Power Systems, Vol. PWRS-2, No. 2, pp. 310-320, May 1987.
[35] V. Venkatasubramanian, H. Schattler, J. Zaborszky, ’Fast time-varying phasor anal-

ysis in the balanced 3-phase large electric power system’, IEEE Transactions on
Automatic Control 40(11): 1975-1982, November 1995.
[36] E.H. Allen, M.D. Ilic, ’Interaction of transmission network and load phasor dy-
namics in electric power systems’, IEEE Transactions on Circuits and Systems I-
Fundamental Theory and Applications 47(11): 1613-1620, November 2000.
[37] P.W. Sauer, B.C. Lesieutre, M.A. PAI, ’Transient algebraic circuits for power system
dynamic modeling’, International Journal of Electrical Power & Energy Systems
15(5): 315-321, October 1993.
[38] P.W. Sauer, M.A. Pai, ’Power System Dynamics and Stability’, Prentice-Hall, 1998.
[39] P.M. Anderson, A.A. Fouad, ’Power System Control and Stability’, The Iowa State
University Press, Volume 1, 1977.
[40] K.R. Padiyar, ’Analysis of Subsynchronous Resonance in Power Systems’, Kluwer

Academic Publishers, 1999.
104 Bibliography
[41] P.C. Krause, ’Analysis of Electric Machinery’, McGraw-Hill, 1987.
[42] P.M. Anderson, B.L. Agrawal, J.E. Van Ness, ’Subsynchronous Resonance in Power
Systems’, IEEE Press, 1990.
[43] W.H. Hayt, J.E. Kemmerly, ’Engineering Circuit Analysis’, McGraw-Hill, 1993.
[44] W.H. Kersting, ’Radial distribution test feeders’, IEEE Power Engineering Society
Winter Meeting, 2001.
[45] http://ewh.ieee.org/soc/pes/dsacom/testfeeders.html.
[46] http://www.abb.com/global/abbzh/abbzh251.nsf!OpenDatabase&db=/global/seitp
/seitp408.nsf&v=9AAC910040&e=us&c=514877510D363EF7C12570180007CEE3
[47] A. Ramirez, A. Semlyen, and R. Iravani, ’Order reduction of the dynamic model of
a linear weakly periodic system-part II: frequency-dependent lines’, IEEE Transac-
tions on Power Systems, Vol. 19, pp. 866-871, May 2004.
[48] S. Carneiro, J.R. Marti, ’Evaluation of corona and line models in electromagnetic
transients simulations’, IEEE Transactions on Power Delivery, Vol. 6, pp.334-342,
January 1991.
[49] A. Morched, B. Gustavsen, and M. Tartibi, ’A universal model for accurate calcula-
tion of electromagnetic transients on overhead lines and underground cables’, IEEE
Transactions on Power Delivery, Vol. 14, pp. 1032-1038, July 1999.
[50] T. Noda, N. Nagaoka, and A. Ametani, ’Phase domain modeling of frequency-

dependent transmission lines by means of an ARMA model’, IEEE Transactions
on Power Delivery, Vol. 11, pp. 401-411, January 1996.
[51] H. V. Nguyen, H. W. Dommel, and J. R. Marti, ’Direct phase-domain modeling

of frequency-dependent overhead transmission lines’, IEEE Transactions on Power
Delivery, Vol. 12, pp. 1335-1342, July 1997.
[52] C. Schauder, H. Mehta, ’Vector analysis and control of advanced static VAr compen-
sators’, Generation, Transmission and Distribution, IEE Proceedings C, Vol. 140, pp.
299-306, July 1993.
Bibliography 105
[53] ’Modeling of power electronics equipment (FACTS) in load flow and stability pro-
grams’, CIGRE TaskForce 38.01.08, August 1999.
[54] C.A. Canizares, ’Power Flow and Transient Stability Models of FACTS Controllers
for Voltage and Angle Stability Studies’, Proceedings of the 2000 IEEE-PES Winter
Meeting, Singapore, January 2000.
[55] H.F. Wang, ’Modelling STATCOM into power systems’, International Conference
on Electric Power Engineering, Power Tech Budapest 99, 1999.
[56] R.Mohan Mathur, R.K. Varma, ’Thyristor-Based FACTS Controllers for Electrical
Transmission Systems’, Wiley-Interscience, 2002.
[57] X.P. Zhang, C. Rehtanz, B. Pal, ’Flexible AC transmission systems: modelling and
control’, Springer-Verlag, 2005.
[58] N. Mohan, T.M. Undeland, and W.P. Robbins, ’Power Electronics, Converters, Ap-
plications and Design (Second Edition)’, John Wiley & Sons, New York, 1995.
[59] Q.H. Wu, J.T. Ma, ’Power system optimal reactive power dispatch using evolution-
ary programming’, IEEE Transactions on Power Systems, Vol. 10, pp. 1243-1249,
August 1995.
[60] K.Y. Lee, X. Bai, and Y.M. Park, ’Optimization method for reactive power plan-
ning by using a modified simple genetic algorithm’, IEEE Transactions on Power
Systems, Vol. 10, pp. 1843-1850, November 1995.
[61] R. Dimeo, K.Y. Lee, ’Boiler-turbine control system design using a genetic algo-
rithm’, IEEE Transactions on Energy Conversion, Vol. 10, pp. 752-759, December
1995.
[62] K.Y. Lee, F.F. Yang, ’Optimal reactive power planning using evolutionary algo-
rithms: a comparative study for evolutionary programming, evolutionary strategy,
genetic algorithm, and linear programming’, IEEE Transactions on Power Systems,
Vol. 13, pp. 101-108, February 1998.
106 Bibliography
[63] E. Diaz-Dorado, J. Cidras, E. Miguez, ’Application of evolutionary algorithms for

the planning of urban distribution networks of medium voltage’, IEEE Transactions
on Power Systems, Vol. 17, pp. 879-884, August 2002.
[64] G. Celli, E. Ghiani, S. Mocci, F. Pilo, ’A multiobjective evolutionary algorithm

for the sizing and siting of distributed generation’, IEEE Transactions on Power
Systems, Vol. 20, pp. 750-757, May 2005.
[65] H. Pohlheim, ’Evolutionre Algorithmen. Verfahren, Operatoren und Hinweise fr die

Praxis’, Springer-Verlag, 1999.
[66] C.M. Hird, H. Leite, N. Jenkins, H. Li, ’Network voltage controller for distributed
generation’, Generation, Transmission and Distribution, IEE Proceedings, Vol. 151,
pp. 150-156, March 2004.
[67] S.N. Liew, G. Strbac, ’Maximising penetration of wind generation in existing distri-
bution networks’, Generation, Transmission and Distribution, IEE Proceedings, Vol.
149, pp. 256-262, May 2002.
[68] N.C. Scott, D.J. Atkinson, J.E. Morrell, ’Use of load control to regulate voltage
on distribution networks with embedded generation’, IEEE Transactions on Power
Systems, Vol. 17, pp. 510-515, May 2002.
[69] N.G. Hingorani, L. Gyugyi, ’Understanding FACTS’, IEEE Press, New York, 1999.
[70] Y.H. Song, A.T. Johns, Eds., ’Flexible AC Transmission Systems (FACTS)’, IEE
Press, London, 1999.
[71] Z. Saad-Saoud, M.L. Lisboa, J.B. Ekanayake, N. Jenkins, G. Strbac, ’Application of

STATCOMs to wind farms’, Generation, Transmission and Distribution, IEE Pro-
ceedings, Vol. 145, pp. 511-516, September 1998.
[72] R. Seydel, ’Practical Bifurcation and Stability Analysis, from Equilibrium to Chaos,
2nd Edition’, Springer Verlag, 1994.
[73] V. Ajjarapu, C. Christy, ’The continuation power flow: a tool for steady state voltage
stability analysis’, IEEE Transactions on Power Systems, Vol. 7, pp. 416 - 423,
February 1992.
Bibliography 107
[74] C.A. Canizares, F.L. Alvarado, ’Point of collapse and continuation methods for large
AC/DC systems’, IEEE Transactions on Power Systems, Vol. 8, pp. 1-8, February
1993.
[75] Hsiao-Dong Chiang, A.J. Flueck, K.S. Shah, N. Balu, ’CPFLOW: a practical tool
for tracing power system steady-state stationary behavior due to load and generation
variations’, IEEE Transactions on Power Systems, Vol. 10, pp. 623-634, May 1995.
[76] X. Zhang, P. Ju, E. Handschin, ’Continuation Three-Phase Power Flow: A Tool

for Voltage Stability Analysis of Unbalanced Three-Phase Power Systems’, IEEE
Transactions on Power Systems, Vol. 20, pp. 1320-1329, August 2005.
Curriculum Vitae
Personal Data:
Born on Oct. 21, 1974, Kunming, P.R. China, Married, Chinese.
Educations:
9/1993 – 7/1997 Bachelor of Science
in the Electrical Engineering Department at Shanghai Jiaotong University,
China.
Besides, the second Bachelor Degree specialized in Mechanical
Engineering.
9/1998 – 7/2001 Master of Science in Engineering

in the Electrical Engineering Department at Kunming University of
Science and Technology, China.
10/2003 – present PhD student

in Lehrstuhls für Energiesysteme und Energiewirtschaft, Dortmund
Universität, Germany.
Work Experience:
8/1997 - 8/1998 Teaching Assistant
Science and Technology, China
8/2001 - 2/2003 Lecturer

Science and Technology, China

Stability of Distribution Systems With A Large Penetration of Distributed Generation

Uploaded by

Copyright:

Available Formats

Stability of Distribution Systems With A Large Penetration of Distributed Generation

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Stability of Distribution Systems With A Large Penetration of Distributed Generation

Uploaded by

Copyright:

Available Formats

Stability of Distribution Systems with a Large Penetration

der Fakultät für Elektrotechnik und Informationstechnik

zur Erlangung des akademischen Grades

Doktor der Ingenieurwissenschaften

Bin Huang M.Sc.(Eng)

Referent: Prof. Dr.-Ing. E. Handschin

Korreferent: Prof. Dr.-Ing. Dr.-Ing. S. Kulig

Tag der mündlichen Prüfung: 9. November 2006

Finally, I wish to express my utmost gratitude to my wife, my parents, my sisters, my

List of the Acronyms iii

2 Concepts, Attributes and Technologies of Distributed Generation and Power

3 Study of the Eigenvalues of the State Space Model of Distribution Networks 14

3.4.1 Properties of Network Eigenvalues . . . . . . . . . . . . . . . . . 27

4 Impacts of DG on the Small Signal Stability of Distribution Systems 39

5 Impacts of DG on the Voltage Stability of Distribution Systems 74

6 Conclusions and Future Works 86

List of the Acronyms

1.1 Motivation and Objectives of this Work

• Possibilities of unsymmetric/unbalanced operations of distribution systems -

• Lack of a general model of DG - Presently there are several technologies for DG

Nevertheless, the question has to be answered: to what penetration level of DG can a

1.2 Outline of the Thesis

Concepts, Attributes and Technologies

2.1 Distributed Generation

2.1.1 Concept and Attributes of Distributed Generation

Nowadays there is a considerable interest in connecting generation to the distribution

Compared to conventional generation, DG usually has characteristics as [15][16]:

• not centrally planned

And it mainly contributes to [15][16]:

• reduction of emission levels (mainly CO2)

2.1.2 Overviews of Distributed Generation Technology

• Microturbines - This new and emerging technology is currently only available

• Industrial combustion turbines - This mature technology ranges from 1 MW to

• Photovoltaics - Commonly known as solar panels, photovoltaic (PV) panels are

2.1.3 Impact of Distributed Generation to Power System Operation

• Protection - The connection of DG changes the magnitude, duration, and direction

quality when it has relatively large capacity or is connected to a weak distribution

2.2 Power System Stability

2.2.1 Overview of Power System Dynamics

Figure 2.1: Time range of power system dynamics [19]

2.2.2 Concepts and Methods for Analyses of Power System Stability

Stability, as a fundamental concept in most engineering disciplines, has been extensively

A more detailed classification according to the different mechanism is helpful to under-

Rotor Angle Stability

Rotor angle stability refers to the ability of synchronous machines of an interconnected

Since the disturbance is assumed to be small, the linearization of system model

2.2.3 Simulation Tools for Analyses of Power System Stability

Study of the Eigenvalues of the State

Main conclusions of this chapter are:

3.1 Modeling of Lines/Cables

3.1.1 Selection of a Common Reference Frame

Figure 3.1: Selection of the common reference frame

The inverse transformation matrix is

3.1.2 Relation of Network Eigenvalues under abc and dq0 Frames

• For λabc2N , assuming it includes complex-conjugate eigenvalues σ1 ± jω1 , after

The discussion consists of two steps.

Then negative sequential sinusoids with frequency of ω1 are assumed to be the

Figure 3.2: System under abc and dq0 frames

λabc2N −18.058 ± j13055 −18.058 ± j13055 −73.708