Manual Of Grid Cal

User Manual:

Open the PDF directly: View PDF PDF.
Page Count: 34

DownloadManual Of Grid Cal
Open PDF In BrowserView PDF
GridCal
Santiago Peñate Vera

G RID C AL
Research oriented power systems software.
Started writing this document in Madrid the 9th of October of 2016

Table of content

1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.1

installation

1.1.1
1.1.2

Windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
Linux / OSX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.2

Run with user interface

1.2.1
1.2.2

Windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
Linux / OSX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.3

Using GridCal as a library

2

Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.1

Circuit and MultiCircuit

2.1.1
2.1.2

Circuit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
MultiCircuit (inherits Circuit) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2

Bus object

2.2.1
2.2.2
2.2.3
2.2.4
2.2.5

Load object . . . . . . . . . . . . .
StaticGenerator object . . . .
Battery object . . . . . . . . . . .
ControlledGenerator object
Shunt object . . . . . . . . . . . .

2.3

Branch object

10

2.4

TransformerType object

11

2.5

Class diagram

11

3

Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.1

Building the admittance matrices

6

6

7

8

9
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

10
10
10
10
10

13

4
3.2

The universal branch model

3.2.1
3.2.2
3.2.3

Yshunt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Yseries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Yf and Yt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.3

The transformer definition from the short circuit test values

3.3.1

Inverse definition of short-circuit values from the Π model . . . . . . . . . . . . . . . . 17

4

Power flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

4.1

Newton-Raphson-Iwamoto

4.1.1

Jacobian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

4.2

Levenberg-Marquardt

20

4.3

DC approximation

22

4.4

Linear AC power flow

23

4.5

Holomorphic Embedding

24

4.5.1
4.5.2
4.5.3
4.5.4

Concepts . . . . . . . . . . . . . .
Fundamentals . . . . . . . . . .
Padè approximation . . . . .
Formulation with PV nodes

4.6

Post power flow (loading and losses)

5

Short circuit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

5.1

3-phase short circuit

6

Graphical User Interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

.
.
.
.

.
.
.
.

14

16

18

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

24
24
27
28

31

32

1. Introduction

GridCal is a research oriented power systems software.
Research oriented? How? Well, it is a fruit of research. It is designed to be modular. As a
researcher I found that the available software (not even talking about commercial options) are hard
to expand or adapt to achieve complex simulations. GridCal is designed to allow you to build and
reuse modules, which eventually will boost your productivity and the possibilities that are at hand.
I have made other projects (even another open source one: fPotencia in C++). I believe that this
project encapsulates my half life of programming experience and the curiosity I have developed for
power systems.
So, I do really hope you enjoy it, and if you have comments, suggestions or just want to
collaborate, do not hesitate to contact.

Cheers,
Santiago

Chapter 1. Introduction

6

1.1

installation
GridCal is designed to work with python 3.5 and onwards, hence retro-compatibility is not guaranteed.

1.1.1

Windows
Open a system console (go to the desktop menu and type cmd, the console should appear). Then
type:
pip install GridCal
This assumes that your system-wide python is a python 3 distribution, hence you must make
sure of this. An easy way to check this, is to open a console, type python and if a python 3 console
opens within the terminal, then it is allright.
For windows systems I have disabled the installation of PyQt (the user interface technology)
because it conflicts with the Qt version provided by Anaconda. This should not be the case, but it
happens. Therefore, under windows, you must install python through Anaconda.

1.1.2

Linux / OSX
On Unix systems the python 2 / 3 issue is non existent since the terminal commands for python 2
and python 3 are different. So, simply go to a system terminal and type:
pip3 install GridCal
This command will install all the dependencies flawlessly unlike under windows.

1.2

Run with user interface
You must have Python 3.5 or higher installed to work with the GUI. From a Python console:
from GridCal.ExecuteGridCal import run
run()

1.2.1

Windows
To run GridCal with GUI directly from a system console, type:
python -c "from GridCal.ExecuteGridCal import run; run()"
You can embed this command in a shortcut file, and have a windows shortcut in the desktop for
GridCal.

1.2.2

Linux / OSX
To run GridCal with GUI directly from a system console, type:
python3 -c "from GridCal.ExecuteGridCal import run; run()"

1.3 Using GridCal as a library

1.3

Using GridCal as a library
You can use the calculation engine directly or from other applications:
from GridCal.Engine.CalculationEngine import *
This will provide access to all the objects in the calculation engine of GridCal.

7

2. Structure

GridCal uses an object oriented approach for all the data and simulation management. However
the object orientation is very inefficient when used in numerical computation, that is why there are
compile() functions that extract the information out of the objects and turn this information into
vectors, matrices and DataFrames in order to have efficient numerical computations. After having
been involved in quite some number-crunching software developments, I have found this approach
to be the best compromise between efficiency and code escalability and maintainability .
The whole idea can be summarized as:
Object oriented structures -> intermediate objects holding arrays -> Numerical modules

2.1

Circuit and MultiCircuit
The concept of circuit should be easy enough to understand. It represents a set of nodes (buses) and
branches (lines, transformers or other impedances)
The MultiCircuit class is the main object in GridCal. It represents a circuit that may contain
islands. It is important to understand that a circuit split in two or more islands cannot be simulated
as is, because the admittance matrix would be singular. The solution to this is to split the circuit
in island-circuits. Therefore MultiCircuit identifies the islands and creates individual Circuit
objects for each of them.
As I said before GridCal uses an object oriented approach for the data management. This
allows to group the data in a smart way. In GridCal I have decided to have only two types of
object directly declared in a Circuit or MultiCircuit object. These are the Bus and the Branch.
The branches connect the buses and the buses contain all the other possible devices like loads,
generators, batteries, etc. This simplifies enormously the management of element when adding,
associating and deleting.

2.1.1

Circuit
• Sbase: Base power to compute the per unit power values (MVA)

2.2 Bus object

9

• branch_original_idx: Array that keeps the index of the branches in the parent MultiCircuit
object.
• branches: List of Branch objects.
• bus_original_idx: Array that keeps the index of the buses in the parent MultiCircuit object.
• buses: List of Bus objects.
• graph: Graph object from the networkx package.
• mc_time_series
• monte_carlo_input
• name: Name of the circuit.
• power_ f low_input
• power_ f low_results
• time_series_input
• time_series_results
2.1.2

MultiCircuit (inherits Circuit)
Since the MultiCircuit inherits the Circuit object, it means that it has by default all the properties of
the Circuit plus the ones that are below.
• branch_dictionary
• branch_names: Array with the name of the branch objects all together.
• bus_dictionary
• bus_names: Array with the name of the bus objects all together.
• circuits: List of the Circuit objects. Each Circuit represents an island.
• has_time_series: Are there time series available? (True/False)

2.2

Bus object
The Bus object is the container of all the possible devices that can be attached to a bus bar or
substation. Such objects can be loads, voltage controlled generators, static generators, batteries,
shunt elements, etc.
• Qmax_sum: Maximum reactive power of this bus (inferred from the devices).
• Qmin_sum: Minimum reactive power of this bus (inferred from the devices).
• V max: Maximum voltage of this bus (p.u.)
• V min: Minimum voltage of this bus (p.u.)
• V nom: Nominal voltage of the bus (kV)
• batteries: List of Battery objects.
• controlled_generators: List of ControlledGenerator objects.
• dispatch_storage: Shall this bus dispatch its storage devices? (True / False)
• graphic_ob ject: Qt graphic object associated to this bus.
• is_enabled: Is this bus enabled for calculation? (True / False)
• is_slack: Is this bus a slack bus? (True / False)
• loads: List of Load objects.
• name: Name of the element.
• shunts: List of Shunt objects.
• static_generators: List of StaticGenerator objects.
• type: Type of the bus.
• x: x coordinate of the bus for representation.
• y: y coordinate of the bus for representation.

Chapter 2. Structure

10
2.2.1

Load object
The load object implements the so-called ZIP model, in which the load can be represented by a
combination of power (P), current(I), and impedance (Z).
• I: Current (complex, in kA)
• I pro f : Pandas DataFrame with the current profile (complex, in kA)
• S: Power (complex in MVA)
• Spro f : Pandas DataFrame with the power profile (complex, in MVA)
• Z: Impedance (complex, in Ohm)
• Z pro f : Pandas DataFrame with the impedance profile (complex, in Ohm)
• name: Name of the load.
The sign convention is: Positive to act as a load, negative to act as a generator.

2.2.2

StaticGenerator object
• S: Power (complex in MVA)
• Spro f Pandas DataFrame with the power profile (complex, in MVA)
• name: Name of the StaticGenerator.

2.2.3

Battery object
•
•
•
•
•
•
•
•
•

2.2.4

ControlledGenerator object
•
•
•
•
•
•
•
•

2.2.5

Enom: Nominal energy capacity (MWh)
P: Active power being dispatched (MW)
Ppro f : Pandas DataFrame with the active power profile (real, in MW)
Qmax: reactive power upper limit (p.u.)
Qmin: reactive power lower limit (p.u.)
Snom: Nominal power (MVA)
V set: Voltage set point (p.u.)
V set_pro f Pandas DataFrame with the voltage set point profile (p.u.)
name: Name of the battery.

P: Active power being dispatched (MW)
Ppro f : Pandas DataFrame with the active power profile (real, in MW)
Qmax: reactive power upper limit (p.u.)
Qmin: reactive power lower limit (p.u.)
Snom: Nominal power (MVA)
V set: Voltage set point (p.u.)
V set_pro f Pandas DataFrame with the voltage set point profile (p.u.)
name: Name of the ControlledGenerator.

Shunt object
• Y : Admittance of the shunt object (in p.u.)
• Y pro f : Pandas DataFrame with the admittance profile (p.u.)
• name: Name of the shunt object.

2.3

Branch object
• angle: Tap angle in radians.
• bus f rom: Bus object to which the branch is connected at the "from" end.
• bust o: Bus object to which the branch is connected at the "to" end.

2.4 TransformerType object
•
•
•
•
•
•
•
•
•

11

is_enabled: Is this branch enabled for calculation?
mtt f : Mean time to failure (h)
mttr: Mean time to repair (h)
name: Name of the branch
rate: Power rating (MVA)
tap_module: tap changer module (value around 1)
y_shunt: Total shunt admittance.
z_series: Total series impedance.
type_ob j: Object of type TransformerType or LineType used to set the current impedance
values.

Implements the model at 3.2.

2.4

TransformerType object
•
•
•
•
•
•
•
•
•

HV _nominal_voltage: High voltage side nominal voltage (kV)
LV _nominal_voltage: Low voltage side nominal voltage (kV)
Nominal_power: Transformer nominal power (MVA)
Copper_losses: Copper losses (kW)
Iron_losses: Iron Losses (kW)
No_load_current: No load current (%)
Short_circuit_voltage: Short circuit voltage (%)
GR_hv1:
GX_hv1:

Implements the model at 3.3.

2.5

Class diagram
Here the API class diagram is sketched. All the classes are included but only the most fundamental
properties and functions od each class are included to keep the diagram simple.

Chapter 2. Structure

12

Battery

<>

MultiCircuit
Bus

+buses: Bus
Inherited from Circuit

+branches: Branch
<>

Inherited from Circuit

+compile()
Overrides

+batteries: Battery
+controlled_generators: ControlledGenerators
+loads: Load
+shunts: Shunt
+static_generators: StaticGenerator

ControlledGenerator

+get_YISV()

Load

returs the total bus admitance, current, power and voltage

Circuit
Shunt

+buses: Bus
+branches: Branch

Branch

+compile()

+apply_to(Ybus:SparseMatrix,Yseries:SparseMatrix,
Yshunt:Vector,i:int,f:int,t:int)

StaticGenerator
Inherits

applies the branch admittance to the relevant matrices and
vectors

Inherits

<>

Line

<>

Transformer

TransformerType

PowerFlowOptions

TimeSeriesInput

VoltageCollapseOptions

PowerFlowInput

<>

VoltageCollapseInput

MonteCarloInput

VoltageCollapseResults

MonteCarloResults

TimeSeriesResults

PowerFlowResults

TimeSeriesResultsAnalysis

<>

<>

<>

PowerFlow

TimeSeries

VoltageCollapse

Figure 2.1: Simplified class diagram of GridCal’s API

MonteCarlo

3. Models

3.1

Building the admittance matrices
This operation occurs in the Compile() function of the Circuit object. This function compiles
many other magnitudes and among them, the following matrices:
• Ybus: Complete admittance matrix.
It is a sparse matrix of size n × n
• Yseries: Admittance matrix of the series elements. It contains no value comming from
shunt elements or the shunt parts of the branch model.
It is a sparse matrix of size n × n
• Yshunt: Admittance vector of the shunt elements and the shunt parts of the branch model.
It is a vector of size n
• Yf: Admittance matrix of the banches with their from bus.
It is a sparse matrix of size m × n
• Yt: Admittance matrix of the banches with their to bus.
It is a sparse matrix of size m × n
Where n is the number of buses and m is the number of branches.
The relation between the admittance matrix and the series and shunt admittance matrices is the
following:
Ybus = Yseries +Yshunt

(3.1)

The algorithmic logic to build the matrices in pseudo code is the following:
n = number of buses in the circuit
m = number of branches in the circuit
For i=0 to n:
bus_shunt_admittance, bus_current, bus_power, bus_voltage = buses[i].get_YISV()
Yshunt[i] = bus_shunt_admittance

Chapter 3. Models

14
end
For i=0 to m:
f = get_bus_inde(branches[i].bus_from)
t = get_bus_inde(branches[i].bus_to)

// the matrices are modified by the branch object itself
branches[i].apply_to(Ybus,Yseries,Yshunt,Yf,Yt,i,f,t)
end

3.2

The universal branch model
The following describes the model that is applied to each type of admittance matrix in the
apply_to() function inside the Branch object seen before.
The model implemented to describe the behavior of the transformers and lines is the π (pi)
model.

Uk
I

t k :1
p

k

Magnetizing
impedance

I
I

bk

Um

q
I

pq

Lekeage
impedance

I

m

bm

Figure 3.1: π model of a branch
To define the π branch model we need to specify the following magnitudes:
• zseries : Magnetizing impedance or simply series impedance. It is given in p.u.
• yshunt : Leakage impedance or simply shunt impedance. It is given in p.u.
• tap_module: Module of the tap changer. It is a magnitude around 1.
• tap_angle: Angle of the tap changer. Angle in radians.
In order to apply the effect of a branch to the admittance matrices, first we compute the complex
tap value.
tap = tap_module · e− j·tap_angle
Then we compose the equivalent series and shunt admittance values of the branch. Both values
are complex.
1
Ys =
zseries
yshunt
Ysh =
2
• zseries : Series impedance of the branch composed by the line resistance and its inductance.
zseries = r + jl

3.2 The universal branch model

15

• yshunt : Shunt admittance of the line composed by the conductance and the susceptance.
yshunt = c + jb


Yf f Yf t
The general branch model is represented by a 2 × 2 matrix. Ybranch =
Yt f Ytt
In this matrix, the elements are the following:
Yf f =

Ys +Ysh
tap · con j(tap)

Y f t = −Ys /con j(tap)
Yt f = −Ys /tap
Ytt = Ys +Ysh
Ybus

The branch admittance values are applied to the complete admittance matrix as follows:
Ybus f , f = Ybus f , f +Y f f
Ybus f ,t = Ybus f ,t +Y f t
Ybust, f = Ybust, f +Yt f
Ybust,t = Ybust,t +Ytt
These formulas assume that there might be something already in Ybus , therefore the right way to
modify these values is to add the own branch values.
3.2.1

Yshunt
Yshunt f = Yshunt f +Ysh
Yshunt t = Yshunt t +

3.2.2

Ysh
tap · con j(tap)

Yseries
Yseries f , f = Yseries f , f +

Ys
tap · con j(tap)

Yseries f ,t = Yseries f ,t +Y f t
Yseriest, f = Yseriest, f +Yt f
Yseriest,t = Yseriest,t +Ys
3.2.3

Yf and Yt
Y f i, f = Y f i, f +Y f f
Y f i,t = Y f i,t +Y f t
Yt i, f = Yt i, f +Yt f
Yt i,t = Yt i,t +Ytt
Here i is the index of the branch in the circuit and f , t are the corresponding bus indices.

Chapter 3. Models

16

3.3

The transformer definition from the short circuit test values
The transformers are modeled as π branches too. In order to get the series impedance and shunt
admittance of the transformer to match the branch model, it is advised to transform the specification sheet values of the device into the desired values. The values to take from the specs sheet are:
•
•
•
•
•
•
•
•

Sn : Nominal power in MVA.
Uhv : Voltage at the high-voltage side in kV.
Ulv : Voltage at the low-voltage side in kV.
Usc : Short circuit voltage in %.
Pcu : Copper losses in kW.
I0 : No load current in %.
GXhv1 : Reactance contribution to the HV side. Value from 0 to 1.
GRhv1 : Resistance contribution to the HV side Value from 0 to 1.

Then, the series and shunt impedances are computed as follows:
2 /S
Nominal impedance HV (Ohm): Znhv = Uhv
n
Nominal impedance LV (Ohm): Znlv = Ulv2 /Sn
Short circuit impedance (p.u.): zsc = Usc /100
Short circuit resistance (p.u.): rsc =
Short circuit reactance (p.u.): xsc =

Pcu /1000
Sn

p

2
z2sc − rsc

HV resistance (p.u.): rcu,hv = rsc · GRhv1
LV resistance (p.u.): rcu,lv = rsc · (1 − GRhv1 )
HV shunt reactance (p.u.): xshv = xsc · GXhv1
LV shunt reactance (p.u.): xslv = xsc · (1 − GXhv1 )
Shunt resistance (p.u.): r f e =

Sn
Pf e /1000

Magnetization impedance (p.u.): zm =
Magnetization reactance (p.u.): xm =

1
I0 /100

r 1
1
2 −
zm

1
r2f e

If the content of the square root is negative, set the magnetization impedance to zero.
The final complex calculated parameters in per unit are:
Magnetizing impedance (or series impedance): zseries = Zm = r f e + j · xm
Leakage impedance (or shunt impedance): Zl = rsc + j · xsc
Shunt admittance: yshunt = 1/Zl

3.3 The transformer definition from the short circuit test values
3.3.1

17

Inverse definition of short-circuit values from the Π model
In GridCal I found the need to find the short circuit values (Pcu , Vsc , r f e , I0) from the branch values
(R, X, G, B). Hence the following formulas:

zsc =

q

R2 + 1/X 2

(3.2)

Vsc = 100 · zsc

(3.3)

Pcu = R · Sn · 1000

(3.4)

zl = 1/(G + jB)

(3.5)

r f e = zl.real

(3.6)

xm = zl.imag

(3.7)

I0 = 100 ·

q
1/r2f e + 1/xm2

(3.8)

4. Power flow

4.1

Newton-Raphson-Iwamoto
The Newton-Raphson method is the standard power flow method tough at schools. GridCal
implements a slight but important modification of this method that turns it into a more robust,
industry-standard algorithm. The Newton Raphson method is the first order Taylor approximation
of the power flow equation. The method implemented in GridCal is the second order approximation,
let’s see how.
The expression to update the voltage solution in the Newton-Raphson algorithm is the following:
Vt+1 = Vt + J−1 (S0 − Scalc )

(4.1)

Where:
• Vt : Voltage vector at the iteration t (current voltage)
• Vt+1 : Voltage vector at the iteration t + 1 (new voltage)
• J: Jacobian matrix.
• S0 : Specified power vector.
• Scalc : Calculated power vector using Vt .
In matrix form we get:

 
  
J11 J12
∆θ
∆P
×
=
J21 J22
∆|V |
∆Q

(4.2)

The formulation implemented in GridCal includes the optimal acceleration parameter µ:
Vt+1 = Vt + µJ−1 (S0 − Scalc )

(4.3)

Here µ is the Iwamoto optimal step size parameter. In 1982 S.Iwamoto and Y.Tamura present a
method [iwamoto1981load] where the Jacobian matrix J is only computed at the beginning, and

4.1 Newton-Raphson-Iwamoto

19

the iteration control parameter µ is computed on every iteration. In GridCal I compute J and µ on
every iteration getting a more robust method on the expense of a greater computational effort.
To compute the parameter µ we must do the following:
Theorem 4.1.1 — Computation of µ.

J’ = Jacobian(Y, dV)
dx = J−1 (S0 − Scalc )
a = S0 − Scalc
b = J × dx
1
c = dx · (J’ × dx)
2
g0 = −a · b
g1 = b · b + 2a · c
g2 = −3b · c
g3 = 2c · c
G(x) = g0 + g1 x + g2 x2 + g3 x3
µ = solve(G(x), x0 = 1)
There will be three solutions to the polynomial G(x). Only the last solution will be real, and
therefore it is the only valid value for µ. The polynomial can be solved numerically using 1 as the
seed.
The matrix J’ is the Jacobian matrix computed using the voltage derivative numerically computed as the voltage increment dV = Vt − Vt−1 (voltage difference between the current and the
previous iteration).

Chapter 4. Power flow

20
4.1.1

Jacobian
The Jacobian matrix is the derivative of the power flow equation for a given voltage set of values.


J11 J12
J=
J21 J22

(4.4)

Where:


• J11 = Re ∂∂ θS [pvpq, pvpq]


• J12 = Re ∂∂|VS | [pvpq, pq]


• J21 = Im ∂∂ θS [pq, pvpq]


• J22 = Im ∂∂|VS | [pqpq]
Where S = V · (I + Ybus × V)∗ .
Here we introduced two complex-valued derivatives:
•
•

∂S
∗
∗
∂ θ = Vdiag · (Ybus ×Vdiag,norm ) + Idiag ·Vdiag,norm
∂S
∗
∂ |V | = 1 j ·Vdiag · (Idiag −Ybus ×Vdiag )

Where:
• Vdiag : Diagonal matrix formed by a voltage solution.
• Vdiag,norm : Diagonal matrix formed by a voltage solution, where every voltage is divided by
its module.
• Idiag : Diagonal matrix formed by the current injections.
• Ybus : And of course, this is the circuit full admittance matrix.
This Jacobian form can be used for other methods.

4.2

Levenberg-Marquardt
The Levenberg-Marquardt iterative method is often used to solve non-linear least squares problems.
in those problems one reduces the calculated error by iteratively solving a system of equations that
provides the increment to add to the solution in order to decrease the solution error. So conceptually
it applies to the power flow problem.
Set the initial values:
• ν =2
• f prev = 109
• ComputeH = true
In every iteration:
• Compute the jacobian matrix if ComputeH is true:
H = Jacobian(Y, V)
• Compute the mismatch in the same order as the jacobian:
Scalc = V(Y · V − I)∗
m = Scalc − S
dz = [Re(m pv ), Re(m pq ), Im(m pq )]

4.2 Levenberg-Marquardt

21

• Compute the auxiliary jacobian transformations:
H1 = H>
H2 = H1 · H
• Compute the first value of λ (only in the first iteration):
λ = 10−3 Max(Diag(H2 ))
• Compute the system Matrix:
A = H2 + λ · Identity
• Compute the linear system right hand side:
rhs = H1 · dz
• Solve the system increment:
dx = Solve(A, rhs)
• Compute the objective function:
f = 0.5 · dz · dz>
• Compute the decision function:
ρ=

f prev −
>

f
0.5 · dx · (λ dx + rhs)

• Update values:
If ρ > 0
– ComputeH = true
– λ = λ · max(1/3, 1 − (2 · ρ − 1)3 )
– ν =2
– Update the voltage solution using dx.
Else
– ComputeH = f alse
– λ = λ ·ν
– ν = ν ·2
• Compute the convergence:
converged = ||dx, ∞|| < tolerance
• f prev = f
As you can see it takes more steps than Newton-Raphson. It is a slower method, but it works
better for ill-conditioned and large grids.

Chapter 4. Power flow

22

4.3

DC approximation
The so called direct current power flow (or just DC power flow) is a convenient oversimplification
of the power flow procedure.
It assumes that in any branch the reactive part of the impedance is much larger than the resistive
part, hence the resistive part is neglected, and that all the voltages modules are the nominal per
unit values. This is, |v| = 1 for load nodes and |v| = vset for the generator nodes, where vset is the
generator set point value.
In order to compute the DC approximation we must perform a transformation. The slack nodes
are removed from the grid, and their influence is maintained by introducing equivalent currents
in all the nodes. The equivalent admittance matrix (Yred ) is obtained by removing the rows and
columns corresponding to the slack nodes. Likewise the removed elements conform the (Yslack )
matrix.

Figure 4.1: Matrix reduction (VD: Slack, PV: Voltage controlled, PQ: Power controlled)

P = real(Sred ) + (−imag(Yslack ) · angle(Vslack ) + real(Ired )) · |Vred |

(4.5)

The equation 4.5 computes the DC power injections as the sum of the different factors mentioned:
1. real(Sred ): Real part of the reduced power injections.
2. imag(Yslack ) · angle(Vslack ) · |vred |: Currents that appear by removing the slack nodes while
keeping their influence, multiplied by the voltage module to obtain power.
3. real(Ired ) · |vred |: Real part of the grid reduced current injections, multiplied by the voltage
module to obtain power.
Once the power injections are computed, the new voltage angles are obtained by:
Vangles = imag(Yred )−1 × P

(4.6)

The new voltage is then:
Vred = |Vred | · e1 j·Vangles

(4.7)

This solution does usually produces a large power mismatch. That is to be expected because
the method is an oversimplification with no iterative convergence criteria, just a straight forward set
of operations.

4.4 Linear AC power flow

4.4

23

Linear AC power flow
Following the formulation presented in [rossoni2016linearized], we obtain a way to solve circuits
in one shot (without iterations) and quite positive results for a linear approximation.

 
 

A11 A12
∆θ
Rhs1
×
=
A21 A22
∆|V |
Rhs2

(4.8)

Where:
• A11 = −Im (Yseries [pqpv, pqpv])
• A12 = Re (Ybus [pqpv, pq])
• A21 = −Im (Yseries [pq, pqpv])
• A22 = −Re (Ybus [pq, pq])
• Rhs1 = P[pqpv]
• Rhs2 = Q[pq]
Here, Ybus is the normal circuit admittance matrix and Yseries is the admittance matrix formed
with only series elements of the π model, this is neglecting all the shunt admittances.
Solving the vector [∆θ + 0, ∆|V | + 1] we get θ for the pq and pv nodes and |V | for the pq nodes.
For equivalence with the paper:
•
•
•
•

−B0 = −Im(Yseries [pqpv, pqpv])
G = Re(Ybus [pqpv, pq])
−G0 = −Im(Yseries [pq, pqpv])
−B = −Re(Ybus [pq, pq])

Chapter 4. Power flow

24

4.5

Holomorphic Embedding
First introduced by Antonio Trias in 2012 [TriasHELM], promises to be a non-divergent power
flow method. Trias originally developed a version with no voltage controlled nodes (PV), in which
the convergence properties are excellent (With this software try to solve any grid without PV nodes
to check this affirmation).
The version programmed in the file HELM.py has been adapted from the master thesis of Muthu
Kumar Subramanian at the Arizona State University (ASU) [subramanian2014application]. This
version includes a formulation of the voltage controlled nodes. My experience indicates that
the introduction of the PV control deteriorates the convergence properties of the holomorphic
embedding method. However, in many cases, it is the best approximation to a solution. especially
when Newton-Raphson does not provide one.
The file HELM.py contains a vectorized version of the algorithm. This means that the execution
in python is much faster than a previous version that uses loops.

4.5.1

Concepts
All the power flow algorithms until the HELM method was introduced were iterative and recursive.
The helm method is iterative but not recursive. A simple way to think of this is that traditional
power flow methods are exploratory, while the HELM method is a planned journey. In theory the
HELM method is superior, but in practice the numerical degeneration makes it less ideal.
The fundamental idea of the recursive algorithms is that given a voltage initial point (1 p.u. at
every node, usually) the algorithm explores the surroundings of the initial point until a suitable
voltage solution is reached or no solution at all is found because the initial point is supposed to be
"far" from the solution.
On the HELM methods, we form a "curve" that departures from a known mathematically
exact solution that is obtained from solving the grid with no power injections. This is possible
because with no power injections, the grid equations become linear and straight forward to solve.
The arriving point of the "curve" is the solution that we want to achieve. That "curve" is best
approximated by a Padè approximation. To compute the Padè approximation we need to compute
the coefficients of the unknown variables, in our case the voltages (and possibly the reactive powers
at the PV nodes).
The HELM formulation consists in the derivation of formulas that enable the calculation of the
coefficients of the series that describes the "curve" from the mathematically know solution to the
unknown solution. Once the coefficients are obtained, the Padè approximation computes the voltage
solution at the "end of the curve", providing the desired voltage solution. The more coefficients
we compute the more exact the solution is (this is true until the numerical precision limit is reached).
All this sounds very strange, but it works ;)
If you want to get familiar with this concept, you should read about the homotopy concept.
In practice the continuation power flow does the same as the HELM algorithm, it takes a known
solution and changes the loading factors until a solution for another state is reached.

4.5.2

Fundamentals
The fundamental equation that defines the power flow problem is:
S = V × (Y × V)∗

(4.9)

4.5 Holomorphic Embedding

25

Most usefully represented like this:
 ∗
S
Y×V =
V

(4.10)

The holomorphic embedding is to insert a "travelling" parameter α, such that for α = 0 we
have an mathematically exact solution of the problem (but not the solution we’re looking for...), and
for α = 1 we have the solution we’re looking for. The other thing to do is to represent the variables
to be computed as McLaurin series. Let’s go step by step.
For α = 0 we say that S = 0, in this way the equation 4.10 becomes linear, and its solution is
mathematically exact. But for that to be useful in our case we need to split the admittance matrix Y
into Yseries and Yshunt . Yshunt is a diagonal matrix, so it can be expressed as a vector instead (no need
for matrix-vector product).
 ∗
S
Yseries × V =
− Yshunt V
V

(4.11)

This is what will allow us to find the zero "state" in the holomorphic series calculation. For
α = 1 we say that S = S, so we don’t know the voltage solution, however we can determine a path
to get there:

Y × V(α) =

αS
V(α)

∗
− αYshunt × V(α) =

αS∗
− αYshunt V(α)
V(α)∗

(4.12)

Wait, what?? did you just made this stuff up??, well so far my reasoning is:
• The voltage V is what I have to convert into a series, and the series depend of α, so it makes
sense to say that V, as it is dependent of α, becomes V(α).
• Regarding the α that multiplies S, the amount of power (αS) is what I vary during the travel
from α = 0 to α = 1, so that is why S has to be accompanied by the traveling parameter α.
• In my opinion the α Yshunt is to provoke the first voltage coefficients to be one. Yseries ×
V[0] = 0, makes V [0] = 1. This is essential for later steps (is a condition to be able to use
Padè).
The series are expressed as McLaurin equations:
∞

V (α) = ∑ Vn α n

(4.13)

n

Theorem 4.5.1 — Holomorphicity check. There’s still something to do. The magnitude

(V(α))∗ has to be converted into (V(α ∗ ))∗ . This is done in order to make the function be
holomorphic. The holomorphicity condition is tested by the Cauchy-Riemann condition, this is
∂ V/∂ α ∗ = 0, let’s check that:

∗

∗

∂ (V(α) ) /∂ α = ∂

∞

∑
n

Vn∗ (α n )∗



∞

/∂ α ∗ = ∑ α nVn∗ (α n−1 )∗
n

(4.14)

Chapter 4. Power flow

26
Which is not zero, obviously. Now with the proposed change:

∗

∗

∞

∗

∂ (V(α )) /∂ α = ∂

∑

V∗n α n



/∂ α ∗ = 0

(4.15)

n

Yay!, now we’re mathematically happy, since this stuff has no effect in practice because our
α is not going to be a complex parameter, but for sake of being correct the equation is now:

Yseries × V(α) =

αS∗
− αYshunt V(α)
V∗ (α ∗ )

(4.16)

The fact that V∗ (α ∗ ) is dividing is problematic. We need to express it as its inverse so it
multiplies instead of divide.
∞
∞
1
= W(α) −→ W(α)V(α) = 1 −→ ∑ Wn α n ∑ Vn α n = 1
V(α)
n=0
n=0

(4.17)

Expanding the series and identifying terms of α we obtain the expression to compute the inverse
voltage series coefficients:

Wn =


1


,


 V0 n

n=0

∑ Wk Vn−k



k=0

 −
,
V0

(4.18)
n>0

Now, the equation 4.16 is:
Yseries × V(α) = αS∗ · W(α)∗ − αYshunt V(α)

(4.19)

Substituting the series by their McLaurin expressions:
∞
n

Yseries × ∑ Vn α = αS
n=0

∗

!∗

∞

∑ Wn α

n

n=0

∞

− αYshunt

∑ Vn α n

(4.20)

n=0

Expanding the series an identifying terms of α we obtain the expression for the voltage
coefficients:

Vn =

0, n = 0
S∗ W∗n−1 −Yshunt Vn−1 ,

n>0

(4.21)

This is the HELM fundamental formula derivation for a grid with no voltage controlled nodes
(no PV nodes). Once a sufficient number of coefficients are obtained, we still need to use the Padè
approximation to get voltage values out of the series.
In the previous formulas, the number of the bus has not been explicitly detailed. All the V
and the W are matrices of dimension n × nbus (number of coefficients by number of buses in the
grid) This structures are depicted in the figure 4.2. For instance Vn is the nth row of the coefficients
structure V.

4.5 Holomorphic Embedding

27
Bus number

n=0

...

n=1

...
...

n=2
n=3

...
.
.
.

.
.
.

.
.
.

.
.
.

.
.
.

.
.
.
...

n=n
Coefficients order

Figure 4.2: Structure of the coefficients
4.5.3

Padè approximation
The equation 4.13 provides us with an expression to obtain the voltage from the coefficients,
knowing that for α = 1 we get the final voltage results. So, why do we need any further operation?,
and what is this Padè thing?
Well, it is true that the equation 4.13 provides an approximation of the voltage by means
of a series (this is similar to a Taylor approximation), but in practice, the approximation might
provide a wrong value for a given number of coefficients. The Padè approximation accelerates
the convergence of any given series, so that you get a more accurate result with less coefficients.
This means that for the same series of voltage coefficients, using the equation 4.13 could give a
completely wrong result, whereas by applying Padè to those coefficients one could obtain a fairly
accurate result.
The Padè approximation is a rational approximation of a function. In our case the function
is V(α), represented by the coefficients structure V. The approximation is valid over a small
domain of the function, in our case the domain is α = [0, 1]. The method requires the function to
be continuous and differentiable for α = 0. Hence the Cauchy-Riemann condition. And yes, our
function meets this condition, we tested it before.
GridCal implements two algorithms that perform the Padè approximation; The Padè canonical
algorithm, and Wynn’s Padè approximation.
Padè approximation algorithm

The canonical Padè algorithm for our problem is described by:

Voltage_value_approximation =

PN (α)
QM (α)

∀α ∈ [0, 1]

(4.22)

Here N = M = n/2, where n is the number of available voltage coefficients, which has to be an
even number to be exactly divisible by 2. P and Q are polynomials which coefficients pi and qi
must be computed. It turns out that if we make the first term of QM (α) be q0 = 1, the function to
be approximated is given by the McLaurin expression (What a happy coincidence!)
PN (α) = p0 + p1 α + p2 α 2 + ... + pN α N

(4.23)

QM (α) = 1 + q1 α + q2 α 2 + ... + qM α M

(4.24)

Chapter 4. Power flow

28

The problem now boils down to find the coefficients qi and pi . This is done by solving two
systems of equations. The first one to find qi which does not depend on pi , and the second one to
get pi which does depend on qi .
First linear system: The only unknowns are the qi coefficients.
qMVN−M+1 + qM−1VN−M+2 + ... + q1VN = 0
qMVN−M+2 + qM−1VN−M+3 + ... + q1VN+1 = 0
...
qMVN + qM−1VN+1 + ... + q1VN+M+1 +VN+M = 0

(4.25)

Second linear System: The only unknowns are the pi coefficients.
V0 − p0 = 0
q1V0 +V1 p1 = 0
q2V0 + q1V1 +V2 − p2 = 0
q3V0 + q2V1 + q1V2 +V3 − p3 = 0
...
qMVN−M + qM−1VN−M+1 + ... + +VN − pN = 0

(4.26)

Once the coefficients are there, you would have defined completely the polynomials PN (α) and
QM (α), and it is only a matter of evaluating the equation 4.22 for α = 1.
This process is done for every column of coefficients V = {V0 ,V1 ,V2 ,V3 , ...,Vn } of the structure
depicted in the figure 4.2. This means that we have to perform a Padè approximation for every
node, using the one columns of the voltage coefficients per Padé approximation.
Wynn’s Padè approximation algorithm

Wynn published a paper in 1969 where he proposed a simple calculation method to obtain the Padè
approximation. This method is based on a table. Weniger in 1989 publishes his thesis where a
faster version of Wynn’s algorithm is provided in Fortran code.
That very Fortran piece of code has been translated into Python and included in GridCal.
One of the advantages of this method over the canonical Padè approximation implementation
is that it can be used for every iteration. In the beginning I thought it would be faster but it turns
out that it is not faster since the amount of computation increases with the number of coefficients,
whereas with the canonical implementation the order of the matrices does not grow dramatically
and it is executed the half of the times.
On top of that my experience shows that the canonical implementation provides a more
consistent convergence.
Anyway, both implementations are there to be used in the code.
4.5.4

Formulation with PV nodes
The section 4.5.2 introduces the canonical HELM algorithm. That algorithm does not include the
formulation of PV nodes. Other articles published on the subject feature PV formulations that
work more or less to some degree. The formulation below is a formulation corrected by myself
from a formulation contained here [liu2017online], which does not work as published, hence the
correction.
Embedding

The following embedding equations are proposed instead of the canonical HELM equations from
section 4.5.2.

4.5 Holomorphic Embedding

29

For Slack nodes:
V (α) = V SP

∀α = 0

(4.27)

For PQ nodes:

 Y × V(α) = 0

∀α = 0

αS
 Y × V(α) = ∗ ∗
V (α )

(4.28)

∀α > 0

For PV nodes:



 Y × V(α) =

S
∀α = 0
V (α ∗ )
S − jQ(α)


∀α > 0
 Y × V(α) =
V∗ (α ∗ )


∗

(4.29)

V (α)V ∗ (α ∗ ) = |V0 |2
∀α = 0
∗
∗
2
SP
2
V (α)V (α ) = |V0 | + (|V | − |V0 |2 ) ∀α > 0

(4.30)

This embedding translates into the following formulation:
Step 1

The formulas are adapted to exemplify a 3-bus system where the bus1 is a slack, the bus 2 is PV
and the bus 3 is PQ. This follows the example of the Appendix A of [liu2017online].
Compute the initial no-load solution (n = 0):
 
  SP 
Vre,1
1
0
0
0
0
0
V [n]re,1
SP 
 0




1
0
0
0
0  V [n]im,1  Vim,1


G21 −B21 G22 −B22 G23 −B23  V [n]re,2   0 

×
=

 B21 G21 B22 G22 B23 G23  V [n]im,2   0 

 
 

G31 −B31 G32 −B32 G33 −B33  V [n]re,3   0 
B31 G31 B32 G32 B33 G33
V [n]im,3
0


∀n = 0

(4.31)

Form the solution vector V[n] you can compute the buses calculated power and then get the
reactive power at the PV nodes to initialize Q[0]:
S = V[0] · (Ybus × V[0])∗

(4.32)

Qi [0] = imag(Si ) ∀i ∈ PV

(4.33)

The initial inverse voltage coefficients W[0] are obtained by:

Wi [0] =

1
Vi [0]

∀i ∈ N

(4.34)

This step is entirely equivalent to find the no load solution using the Z-Matrix reduction.

Chapter 4. Power flow

30
Step 2

Construct the system of equations to solve the coefficients of order greater than zero (n > 0). Note
that the matrix is the same as constructed for the previous step, but adding a column and a row for
each PV node to account for the reactive power coefficients. In our 3-bus example, there is only
one PV node, so we add only one column and one row.


 
 

1
0
0
0
0
0
0
V [n]re,1
0
 0

 

1
0
0
0
0
0 

 V [n]im,1   0 
G21 −B21 G22




−B22 G23 −B23 W [0]im  V [n]re,2   f 2re 


 B21 G21

 

B22
G22
B23 G23 W [0]re 

 × V [n]im,2  =  f 2im 
G31 −B31 G32

 

−B32 G33 −B33
0 

 V [n]re,3   f 1re 
 B31 G31




B32
G32
B33 G33
0
V [n]im,3
f 1im 
0
0
V [0]re V [0]im 0
0
0
Q2 [n]
ε[n]

∀n > 0 (4.35)

Where:
f 1 = Si∗ ·Wi∗ [n − 1] ∀i ∈ PQ

(4.36)

f 2 = Pi ·Wi∗ [n − 1] + conv(n, Qi ,Wi∗ ) ∀i ∈ PV

(4.37)

 1
1
|ViS P|2 − |Vi [0]|2 − conv(n,Vi ,Vi∗ ) ∀i ∈ PV, n > 0
2
2
The convolution conv is defined as:
ε[n] = δn1 ·

(4.38)

n−1

conv(n, A, B) =

∑ A[m] · B[n − m]

(4.39)

m=0

The system matrix (Asys ) is the same for all the orders of n > 0, therefore we only build it once,
and we factorize it to solve the subsequent coefficients.
After the voltages V[n] and the reactive power at the PV nodes Q[n] is obtained solving the
linear system (eq 4.35), we must solve the inverse voltage coefficients of order n for all the buses:
n

− ∑ Wi [m] ·Vi [n − m]
Wi [n] =

m=0

Vi [0]

∀i ∈ N, n > 0

(4.40)

Step 3

Repeat step 2 until a sufficiently low error is achieved or a maximum number of iterations (coefficients).
The error is computed by comparing the calculated power S (eq 4.32) with the specified power
injections SSP :
mismatch = S − SSP

(4.41)

error = |mismatch|∞ = max(abs(mismatch))

(4.42)

4.6 Post power flow (loading and losses)

4.6

31

Post power flow (loading and losses)
As we have seen, the power flow routines compute the voltage at every bus of an island grid.
However we do not get from those routines the "power flow" values, this is the power that flows
through the branches of the grid. In this section I show how the post power flow values are computed.
First we compute the branches per unit currents:

If = Yf ×V

(4.43)

It = Yt × V

(4.44)

These are matrix-vector multiplications. The result is the per unit currents flowing through a
branch seen from the from bus or from the to bus.
Then we compute the power values:
S f = V f · I∗f

(4.45)

St = Vt · It∗

(4.46)

These are element-wise multiplications, resulting in the per unit power flowing through a branch
seen from the from bus or from the to bus.
Now we can compute the losses in MVA as:

losses = |S f − St | · Sbase

(4.47)

And also the branches loading in per unit as:
loading =

max(|S f |, |St |) · Sbase
rate

The variables are:
• Y f , Yt : From and To bus-branch admittance matrices, see section 3.1.
• I f : Array of currents at the from buses in p.u.
• It : Array of currents at the to buses in p.u.
• S f : Array of powers at the from buses in p.u.
• St : Array of powers at the to buses in p.u.
• V f : Array of voltages at the from buses in p.u.
• Vt : Array of voltages at the to buses in p.u.
• rate: Array of branch ratings in MVA.

(4.48)

5. Short circuit

5.1

3-phase short circuit
First, declare an array of zeros of size equal to the number of nodes in the circuit.
I = {0, 0, 0, 0, ..., 0}
Then, compute the short circuit current at the selected bus i and assign that value in the ith
position of the array I.
Ii = −

V pre− f ailure,i
Zi,i + z f

(5.1)

Then, compute the voltage increment for all the circuit nodes as:
∆V = Z × I

(5.2)

Finally, define the voltage at all the nodes as:

V post− f ailure = V pre− f ailure + ∆ V

(5.3)

Magnitudes:
• I: Array of fault currents at the system nodes.
• V pre− f ailure : Array of system voltages prior to the failure. This is obtained from the power
flow study.
• z f : Impedance of the failure itself. This is a given value, although you can set it to zero if
you don’t know.
• Z: system impedance matrix. Obtained as the inverse of the complete system admittance
matrix.

6. Graphical User Interface

Figure 6.1: GridCal user interface representing the IEEE 39 bus grid split in three islands.

Bibliography
Source Exif Data: 
File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.5
Linearized                      : No
Page Count                      : 34
Page Mode                       : UseOutlines
Author                          : Author
Title                           : Title
Subject                         : 
Creator                         : LaTeX with hyperref package
Producer                        : pdfTeX-1.40.18
Create Date                     : 2018:10:12 11:43:36+02:00
Modify Date                     : 2018:10:12 11:43:36+02:00
Trapped                         : False
PTEX Fullbanner                 : This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) kpathsea version 6.2.3
EXIF Metadata provided by EXIF.tools

Navigation menu