Developing a harmonic power flow software in distributed generation systems

(1)

Degree project in

CÉDRIC BUREAU

Stockholm, Sweden 2012

XR-EE-ES 2012:011 Electric Power Systems

Second Level,

(2)

Accessibilité : Page 1 sur 58 EDF

Master thesis report

Developing a harmonic power flow software in distributed generation

systems

Student: Cédric BUREAU Supervisor: Lennart SÖDER

EDF R&D, Département Mesure et Système d’Information des Réseaux Electriques 1, avenue du Général de Gaulle

92141 CLAMART

(3)

Abstract

The main topic of this thesis is harmonic power flow and its use in a simulation software that I have developped. The idea of the software is to combine distribution grids’

description, non-linear load models and power flow methods.

Nowadays, power electronics is more and more present in electric devices in distributed generation systems. Those power electronics systems can emit or absorb harmonics that can damage the devices in the grid. Thus, it is important to be able to estimate harmonic behaviour in the grid in order to be able to prevent the possible problems that could occur.

The main contribution of this internship is the precise expression of the needs and goals of the software, and an implementation of its structure. In this thesis, It is explained how the grid’s components and non-linear devices are modelled in the software in order to be able to represent the distribution system. There is also a study the possible input of this software and create a symbolic representation of the grid that is helpful when it comes to load flow calculation. Then, the different load flow and harmonic load flow algorithms that are presented in the literature are analysed and compared them together in order to determine the methods that should be implemented in the future software.

Two of the implemented fundamental load flows with a single-phase system are tested.

Thus, it also validates the input reading and the grid representation construction.

The software developped is a first implementation of a more global software that will require further studies. Indeed, the developpement stage will be done by external

contractors or computer science specialists, that will insist on parallelization of algorithms

and software optimization, in order to have a software as efficient and fast as possible.

(4)

Acknowledgements

First of all, I would like to express my gratitude to my supervisor Professor Lennart Söder, for his advice and support during the project.

I would like to thank the members of the group Power Quality of EDF R&D and especially

Vincent Diemunsch and Christophe Santander for their support and encouragement and for

providing stimulating atmosphere during my internship.

(5)

1. Introduction

1.1. Presentation of Power Quality Group in MIRE at EDF

1.1.1. Presentation of the company

EDF (« Electricité De France ») is one of the world’s largest utility company and is present in every parts of the electricity activities, from production to distribution. Since the promulgation of the law in 2004 about the deregulation of the electricity market, two major activities have become subsidiaries of EDF: the electricity power transmission managed by RTE and the distribution managed by ErDF.

The group exists mainly in Europe and South America and tries to develop its activities in Asia and in the USA.

Some relevant statistics:

- around 160000 employees - 37 million of customers - 2009 Turnover: €63.34 billion - 630.4 TWh produced in 2010

1.1.2. Presentation of Power Quality Group

There are more than 2 000 researchers in the Research and Development department in EDF.

The main mission of the department “Measurement and Information Systems of Electrical Networks” is to analyze the behavior of the transmission system, the distribution system, or even private electrical systems. Employees try to increase the global efficiency of those networks. Its annual budget is 13 M€.

The Power Quality Group’s mission is to develop methods and tools that will help to control that the electricity sold by EDF respect the quality agreement and the needs of the customers. Some of its activities are:

- representing EDF in the national and international standardization institution - offering an expertise about electricity quality to other groups in EDF or ErDF - developing or testing new instrumentation solution

- developing and maintaining the computer tools that helps to realize studies about electricity quality

1.2. Background

Power electronics is more and more present in electric domestic and industrial devices, in distributed generation systems and even in electric cars chargers. Those power electronics systems have a tendency to absorb or emit distorted currents from a 50 Hz sinusoid. Thus, harmonic voltages appear in the grid and they can damage the systems.

In order to prevent the problems, it can be interesting to be able to simulate existing network behavior when some non-linear buses are connected and observe where harmonics voltages appear.

The combination of non-linear loads’ models and harmonic power flow software could give interesting information in order to improve systems stability and safety.

In the Power Quality Group, non linear loads frequency models are built. Furthermore, a module that realizes time domain simulation of non-linear loads from those models is developed. Thus, there is a need for a simulation software that could simulate existing grids and realize a harmonic power flow by communicating with that module in order to get an accurate response of the loads.

This internship’s goal is to write the specifications, develop such a software and compare the results

(9)

Accessibilité : Page 8 sur 58 EDF with existing simulation methods.

1.3. Plan of the thesis

In the following part of this thesis, we will globally describe the software we are developing by presenting its purpose and how it operates. Then, we will explain the contribution of this software as regards the existing software used by EDF.

In the third part, we will describe the components’ library by introducing some components models.

In the fourth part, we will describe the possible input format and the output format that would be interesting to develop.

In the fifth part, we will explain how the network’s description will be internally represented and how it will be used by the calculation module.

In the following part, we will describe the existing algorithm that can be used in order to realize fundamental and harmonic power flow. We will try to select some of them to be implemented in the software.

In the seventh part, we will briefly explain the software’s architecture by presenting the uses cases and the main packages that constitute the software.

In the last part, we will discuss the results by comparing our software with the existing ones.

2. Software description

2.1. Global description

The goal of this software is to read a real electrical network description and to make determinist frequency calculations in order to obtain harmonic levels. It must take into account three-phase networks, distributed generation and non-linear loads.

This software shall also be able to communicate with an external module that is developed by the department in order to realize time domain simulation of non-linear loads. For this study, we will consider that the network is linear. Furthermore, I will not implement the connection with the external module. Thus, we will only use known models of non-linear loads.

We shall also be able to be connected with OpenTURNS, a software that models uncertainty of parameters, in order to realize multiple simulations of the same network, modifying at each iteration the value of a few parameters.

(10)

Figure 1 - Software description

In this first implementation of the software, I will not implement the connection with OpenTURNS. We will modify ourselves the network if we want to. The user will be able to choose what kind of load flow he wants to use and the voltages and currents he wants to observe.

Even if all the functionalities are not implemented in the first version, we will take care of developing the software in the most modular way in order to easily integrate future modules.

2.2. Definition of basic elements

2.2.1. Node and Bus definition

A node represents a physical point in the system.

A bus represents a group of nodes to which a component is connected.

For instance we have for a three-phase component:

Figure 2 - Three-phase component example

2.2.2. Harmonic definition

A harmonic is a component frequency of the signal that is a multiple of the fundamental frequency

Hz f

₀

= 50

.

Thus, we call the

h

^th harmonic the signal that has a frequency

f = h . f

₀. It is defined by the equation:

{ } ⁿ

i ∈ 1 ..

∀

^,

^V

ⁱ^h

⁼ ^V

ⁱ^h

^. ^exp ( ) ^j ^. ^θ

ⁱ^h

The voltage angle

θ

refers to the reference voltage which is taken at the source of the studied system, which is usually the point before the transformer of the substation.

In this study, the last harmonic we will take into account is the

L

^th harmonic.

2.3. Contribution of this software

2.3.1. Existing software

2.3.1.1. HARMONIQUE

The software HARMONIQUE was developed by EDF and is used to calculate frequency simulation. It can calculate voltages and currents on LV and MV electrical networks.

This software has many issues:

(11)

Accessibilité : Page 10 sur 58 EDF - It can only simulate single-phase grids.

- The calculation is independent for every harmonics.

- The current sources used to model non-linear loads are calculated without taking into account the voltage at the considered node.

- The grid size is limited to 300 nodes.

- There is no load flow calculation.

2.3.1.2. ExperTEC

This software was also developed by EDF in the department I am doing this internship. It can realize permanent and transcient analysis, and calculate harmonic voltages values.

The main advantages of this software are:

- It can simulate three phase unbalanced systems.

- Components are very detailed and their parameters can be easily changed by a user.

- It can realize simulation on networks that contains up to 5000 nodes.

However, this software has many issues:

- It does not realize a real load flow.

- It has been developed by one person that left the department.

2.3.2. Advantages of the software

The main advantage of our new software is that it will realize both fundamental and harmonic load flow. Different methods will be implemented in order to compare algorithms and determine the most efficient, by taking into account the accuracy of the results and the execution time.

Besides, we will be able to connect this software with an external module that realizes temporal simulation of a non-linear load. Thus we could use this module at each iteration of the harmonic load flow to get an accurate value of the non-linear load behavior, corresponding to a value of harmonic voltages at that load.

Furthermore, we will also be able in the future to use OpenTURNS in order to realize parametric simulations. The software and the intermediary file format will be designed in order to ease the components’ or network’s modifications done by OpenTURNS.

Finally, this software will be able to read files that describes existing distribution network. Thus, we could observe the impact of a non-linear load connection on a specific real grid.

3. Library of components - frequency models

3.1. Library presentation

In order to be able to realize a harmonic load flow, we need to have a frequency model for every network components and ending components called equipments.

We choose to calculate the models from the physical or basic electrical data usually given by the components manufacturers. The input data format shall give those data.

(12)

Accessibilité : Page 11 sur 58 EDF The idea is that we should be able to calculate the

Y

_busmatrix from the manufacturers’ data.

The components we take into account can be single-phase or three-phase electrical components.

Furthermore, we describe the basic elements Z that can be used to model losses in a component.

Figure 3 - components library

In the following section, we will describe the complex calculation of one components

Y

_busmatrix for each frequency, which take into account impact of the frequency on the component behavior.

The components

Y

_busmatrix describes the component behaviour. For example, for the component

described in Figure 2, the matrix is described by the equation

 





 





=

 





 





6 5 4 3 2 1

.

V V V V V V

Y

I I I I I I

bus .

In order to write the equations, I used the work made previously in order to develop the software ExperTEC. The development of those models has been done and described in an internal report concerning ExperTEC (cf. [17], [18], [19]).

3.2. Matrices for each complex component

In order to model the components, we have to be able to get a component admittance matrix for each

(13)

Accessibilité : Page 12 sur 58 EDF component of the network and for each frequency. Thus, we will consider in the following part that the study is done at a given frequency

f = h . f

₀^.

3.2.1. Model utilization

Here we will explain how we can use simple models of single-phase cable, line in order to obtain three phases models of those components. For transformer, we can either use three single phase model to calculate the three phase component admittance matrix or use directly a three phase model based on the electromagnetic equations and then take into account the impact of windings coupling.

3.2.1.1. Line and cable

Figure 4 – Line and cable model

The p and s indexes represent respectively the primary and secondary parts of the component.

In order to analyze this element, we will study separately its serial part and its parallel part. First, we model the serial part of the line.

With this model, the impedance matrix at frequency f of the serial part is, for the abc components:

 







 







 







 







=

 







 







−

=

 







 







∆

c b a

p m m

m p m

m m p

sc pc

sb pb

sa pa

c b a

I I I

Z Z Z

V V

V V V

.

The entries of the matrix depend on the frequency f.

Then, we use the symmetrical components and the component impedance matrix at frequency f becomes:

 







 







 







 







−

− +

=

 







 







∆

2 1 0

0 0

2 I I I

Z Z Z Z Z Z

V V V

m p m p m p

Then, we model the parallel part of the line (the shunt impedance).

For the abc parameters, the component admittance matrix at frequency f of the parallel part becomes:

(14)



 







 





 







 







+

−

− +

−

− +

=



 







 





pc pb pa

T T

T

pc pb pa

V V V

Y Y

I I I

. . 2 .

2 .

2

ϕϕ ϕ

ϕϕ ϕϕ

ϕ ϕϕ

ϕϕ ϕϕ

ϕϕ ϕ

.

Then, we use the symmetrical components and the admittance matrix at frequency f becomes:

 







 







 







 







+ +

=

 







 







2 1 0

. . 3 0

0 0 .

3 0

0 0

p p p

T T

T

p p p

V V V

Y Y

Y

I I I

ϕϕ ϕ

ϕ

.

Thus, we can consider for the symmetrical components that a three phase line is actually three independent single-phase lines. The parameters of those lines can be calculated from the parameters of the three-phase line.

We can indeed calculate for each symmetrical component the admittance matrix of the corresponding line, at frequency f:

 

 



 



 





+

−

= +

 

 





s p s p s

s s

p

s p

V V Y Y Y

Y Y

Y I

I .

Thus we can determine the component admittance matrix for all the symmetrical components, at frequency f, by combining the three admittance matrices of the three symmetrical components:

 





 





+

−

+

−

+

−

− +

=

s p s

s s

p

s s

p

s s

p

Y Y Y

Y Y

Y

Y Y

Y

Y Y

Y

2 2 2

1 1 1

0 0 0

2 2

2

1 1

1

0 0

0

012

0 0

By taking into account the length of the line, we have:

Y

_T_{_}₀₁₂.

The component admittance matrix for the abc components for the differential values of the voltages at frequency f is then obtained:

1 012

_

0 . 0 0 .

0

⁻

 



 



 



 

= 

S S T

S S

abc

T

Y T T

Y T

Where

T

_S is Fortescue’s Matrix:

 







 







=

2 2

1 1

1 1 1

α α

T

S .

3.2.1.2. Three-phase transformer First Method

In this part, we will explain how we can use a model of single-phase transformer that is different for zero sequence and positive, negative sequence, in order to calculate the component admittance matrix of a three-phase transformer.

The phases of the transformer are called a, b and c and we have:

(15)

[ ]   





 







=

c p b p a p abc p

I I I

I

,

[ ]

 







 







=

c s b s a s abc s

I I I

I

,

[ ]

 







 







=

c p b p a p abc p

V V V

V

and

[ ]

 







 







=

c s

b s

a s abc s

V V V V

The component admittance matrix of the transformer at frequency f is defined by:

 

 



= 

 

 





abc s

abc p abc

s abc p

V Y V I

I .

with





 



= 

SS SP

PS PP

Y Y

Y Y Y

Then, we use the symmetrical component:

[ ] [ ] ^V

^p^abc

⁼ ^T

^S

^. ^V

^p⁰¹² ^and

[ ] [ ] ^I

^abc^p

⁼ ^T

^S

^. ^I

⁰¹²^p ^with

[ ]   





 







=

2 1 0 012

p p p

p

I I I

I

,

[ ]

 







 







=

2 1 0 012

s s s

s

I I I

I

,

[ ]

 







 







=

2 1 0 012

p p p

p

V V V

V

and

[ ]

 







 







=

2 1 0 012

s s s

s

V V V V

The component admittance matrix at frequency f becomes:

 





 





 =



 



= 

⁻₋ ⁻₋

2 2

1 1

0 0

2 2

1 1

0 0

1 1

012

0 0 0

0 0

0 0 0

0 0 0 0

0 0 0

0 0

0 0 0

0 . . .

.

. . .

.

ss sp

ps pp

S SS S S SP S

S PS S S PP S

Y Y

T Y T T Y T

T Y T T Y

Y T

.

We obtain three independent modes. We can determine for each of them a value of the admittance matrix by using a model of single-phase transformer. Entries of matrices depend on the frequency f.

Thus, by calculating the value of the zero, positive and negative sequence, we can obtain the values of

Y

_PP,

Y

_SP,

Y

_PSand

Y

_SSand then the matrix

Y

for the abc components, at frequency f.

Second Method

From the transformer electromagnetic equations, we use a model described in an internal report within EDF R&D (cf. [17]) for each kind of flux: free flux, independent flux and forced flux.

We simplify those models with the following hypothesis: equality of reluctances. It means that the transformer has a symmetric structure. Thus, as it is described in [17], we find a 6x6 matrix

Z

_abc that is the differential impedance matrix in the abc components system, that describes the relation between

currents and potential differences, at frequency f:

 



 







 



 







=

 



 







 



 







sc sb sa pc pb pa

abc

sc sb sa pc pb pa

I I I I I I

Z

V V V V V V

.

The entries of

Z

_abc depend on the frequency f.

By inverting it, we find the corresponding admittance matrix

Y

_abc.

(16)

Figure 5 – Differential and nodal representations of transformer

In order to take into account the impact of windings coupling on the component admittance matrix, we have to calculate the nodal admittance matrix

Y

_Nodal_{_}_abc(12x12) at frequency f defined by

 



 







 



 







=

 



 







 



 







11 10 1 0

_

11 10 1 0

...

. ...

V V V V

Y

I I I I

abc

Nodal , by doing the following calculation:

t

abc abc

Nodal

Y

 



 





 



 





−

 



 





 



 





−

=

1 1 1 1 1 1

0

0 1 1 1 1 1 1

.

1 1 1 1 1 1

0

0 1 1 1 1 1 1

_

Then, we calculate a matrix

C

that represents the impact of windings coupling on this admittance matrix. For instance, for a Wye-Delta-1 transformer, this coupling matrix will be:

(17)

 



 





 



 





=

0 0 0 0 0 0

0 1 1 0 0 0

0 0 0 1 1 0 0

1 0 0 0 0 1 1 0 1 0 1 0

0 1 0 0 0 0

0 0

0 0 1 0 0

0 0 0 0 0 1

C

Finally, we get an 8x8 component admittance matrix

Y

that fully describes the transformer behavior at

a given frequency f:



 



 









 



 







=



 



 









 



 







s c

s b

s a

s neutral

p c

p b

p a

p neutral

s c

s b

s a

s neutral

p c

p b

p a

p neutral

V V V V

Y

I I I I

, , ,

, , , ,

,

, , ,

, , , ,

,

.

The calculation is made using the following formula:

Y = C . Y

_Nodal_{_}_abc

C

^t. It is this matrix that will be returned when one asks its admittance at a given frequency to a component.

3.2.2. Example of models calculation

3.2.2.1. Single-phase transformer

This study is done at a given frequency

π ω

.

= 2

f

.

Manufacturer’s data:

V

p : primary rms voltage in V

V

s : secondary rms voltage in V

S

n : apparent power in VA

v

sc : short circuit voltage in V

I

0 : no-load current in A

R

p,

R

_s : winding resistance in

Ω / m

Model:

(18)

Figure 6 – Single-phase transformer model Model’s elements calculation:

First we have to calculate:

The turns ratio:

s p

V n = V

primary leakage inductance:

ω . 2

11

100

sc

Z

np

l = v

with

n p

np

S

Z V

2

=

secondary leakage inductance:

ω . 2

22

100

ns

sc

Z

l = v

^with

n s

ns

S

Z V

2

=

Thus, we can know the model elements:

.I

0

L

_p

V

^p

= ω

L

p

k = 1 − l

¹¹

k L

_s

l

= − 1

22

s p

L L k M = . .

Ybus matrix, at frequency

π ω

.

= 2

f

:

 

 



= 

p m

m s

y y

y

Y y

with:

( ) ( )

 





 







− + +

+

= −

− + +

+

= +

− + +

+

= +

2 2

. .

.

. .

.

. .

.

. .

ω

ω ω ω

ω

ω ω

s p p

s s p s

p m

s p p

s s p s

p

p p p

s p p

s s p s

p

s s s

L L M L

R L R j R R

M y j

L L M L

R L R j R R

L j y R

L L M L

R L R j R R

L

j

y R

(19)

Accessibilité : Page 18 sur 58 EDF 3.2.2.2. Single-phase autotransformer

This study is done at a given frequency

π ω

.

= 2

f

.

Data: the given data are the same as transformer’s one.

Model:

Figure 7 – autotransformer model Model’s elements calculation

Calculations of

M

,

k

,

L

_pand is the same as those we explained for the transformer. However, we have to calculate

R

_poand

L

_po:

s p

po

R R

R = −

L

po is calculated by solving the following equation:

^L

^po²

⁻ ² ^. ^L

^po

( ^L

^p

⁺ ( ² ^k

²

⁻ ¹ ) ^L

^s

) ( ⁺ ^L

^p

⁻ ^L

^s

)

²

⁼ ⁰

Ybus matrix, :

Increasing voltage auto transformer:

( ) ( )

(

^p

)

^p ^p ^p ^so ^s ^so

(

^p ^so

)

p

p p

so p so

p

L L M R

R j R

L R j R M

L j R

M L j R M L

L j R Y R

. .

. . . . 2 . 1

2

−

+

 +



 





+ +

−

+

−

− +

+ +

= +

ω ω ω

ω

ω ω

Decreasing voltage autotransformer:

( )

(

^s

)

^s ^po

(

^po ^s

)

^s ^po ^p ^po

(

^s ^po

)

s

s s

s p

L L M R

R j R

M R L L j R R M L j R

M L j R L

j Z R

. .

. . . . 2 . 1

2  + +

² ²

−





 



+ + +

+ +

−

+

−

= +

ω ω ω

ω

ω ω

3.2.2.3. Three phase cable

In order to model, the three phase cable, I use model that have been developed in EDF R&D. Those models are confidential and I cannot detail them in the report.

(20)

Accessibilité : Page 19 sur 58 EDF With this model, we can find positive, negative and zero sequence a numerical value of admittance matrix, at frequency f.

 



 

 





 



+

−

= +

 



 



s p s p s

s s

p

s p

V V Y Y Y

Y Y

Y I

I .

Thus, we can have the component admittance matrix for the three-phase component at frequency f as it is explained in 3.2.1.1:

1

2 2 2

1 1 1

0 0 0

2 2

2

1 1

1

0 0

0

0 . 0

0 0

0 .

0

⁻

 

 





 





 





+

−

+

−

+

−

− +

 

 



= 

S S

s p s

s s

p

s s

p

s s

p

S S

abc

T

Y Y Y

Y Y

Y

Y Y

Y

Y Y

Y

T Y T

3.2.2.4. Three-phase line

In order to model, the three-phase line, I use model that have been developed in EDF R&D. Those models are confidential and I cannot detail them in the report.

With this model, we can find for the positive, negative and zero sequence a numerical value of admittance matrix.

 



 

 





 



+

−

= +

 



 



s p s p s

s s

p

s p

V V Y Y Y

Y Y

Y I

I .

Thus, we can have the component admittance matrix for the three-phase component at the frequency f as it explained in 3.2.1.1.

3.2.2.5. Impedance

This component is a fixed impedance and the frequency does not modify its value.

Data:

R

: resistance in

Ω X

: reactance in

Ω

Ybus matrix:

- If the impedance is connected between two different nodes of the network, the component

admittance matrix is:



 









 







+

− +

= +

X j R X j R

X j R X

j Y R

. 1 .

1 . 1 .

1

.

- If the impedance is connected between a node of the network and the ground, the component admittance matrix that is used is:







 



= +

X j Y R

.

1

.

(21)

3.2.3. Linear Bus

In this study, we only consider single-phase linear loads.

Those components are not included in the system admittance matrices, for every harmonics.

At fundamental frequency, linear buses are defined by the two parameters that are known. Thus, a linear bus can be a PQ bus, a PV bus or a

^V θ

bus which is called the slack bus.

Even if those buses do not produce harmonic currents, they have a passive behavior as regards harmonics that go through the network. This behavior can be modeled by the following impedance (cf.

[3]), where the upper part is an impedance and the lower part is a reactance.

Figure 8 – Linear Bus Harmonic Model

With

P

¹: fundamental active power

Q

¹: fundamental reactive power

V

¹: fundamental voltage

h

: harmonic rank

This model is equivalent to a current injection, which is the sum of the current injection in the reactance and the current injection in the impedance:

( ) ( )

¹ ²

1 1 2

1

. . . .

V h

Q V j V

P I V

h h

h

= −

3.2.4. Non linear Bus

Those components are not included in the system admittance matrix.

The fundamental behavior of those buses is the same behavior as PQ bus. It means that during a fundamental load flow calculation, we will consider that the NLD has a known active and reactive fundamental power.

For harmonic calculations, we will consider that the non-linear loads are equivalent harmonic current injectors (cf. [2], [23]). Models that describe those loads are obtained by different ways:

- Some models that describe very specific loads have been developed by a partner university for EDF (cf. [21]). The format of injected current according to those models is given by:

(

ⁱ ⁱ ⁱ^L ^eq

)

k k

i

g V V V Z

I =

¹

,

²

,.., ,

where

L

: last harmonic

k

: harmonic number

i

: node number

( )

1 1 2

P V

( )

1 1 2

.

Q

V

h

Developing a harmonic power flow software in distributed generation systems

Degree project in

CÉDRIC BUREAU

Stockholm, Sweden 2012

XR-EE-ES 2012:011 Electric Power Systems

Second Level,

Master thesis report