Power flow solver for power systems analysis

A function that solves the power flow equations of an electrical power network

Two files are included in this submission:

power_flow_solver.m – is the solver itself (matlab function)

test_flow_solver.m – Example file. This script defines a simple power network and uses 'power_flow_solver' function to compute the corresponding power flow.

Just download the two files to your Matlab directory and it should work fine.

Full documentation may be accessed by typing:
'HELP power_flow_solver' at the matlab command prompt.

Written by Dr. Yoash Levron, October 2012.
I work at the Technion University (Haifa, Israel).
my E-mail: [e-mail]
You can contact me regarding questions.
Use of this function is free of charge, however references will be appreciated in case it supports published research.

Power-flow solver are extensively used in power system analysis. This function provides a basic power-flow solver, based on the Gauss-Seidel method.
The function may solve a power network of arbitrary size and complexity. The solution is fast and robust.
Convergence is guaranteed if the network does have a solution.
However, this function is limited to the analysis of P-Q buses only (representing either generators or loads). P-V buses cannot be modeled. In addition, voltage and current limits are ignored. Therefore, this function may be used
in the analysis of small power networks, or as an instructive tool for students, but is probably impractical in real-world
power networks.

The network is composed of power lines (admittances) and buses.
buses are indexed 1..N. A power source and/or a shunt admittance may be connected to each bus. The direction of power is always defined towards the network. Generators are defined with positive real power, and loads with negative real power.
All buses are considered P-Q buses, with one exception: Bus 1 is always the 'slack' bus, represented as voltage source having a preset phase (this phase serves as a reference for all other phases). The usual choice is E1 = 1.0*exp(j*0), that is, a slack bus with unity voltage and zero phase.

function [Ebus, Ibus, Imat, iter] = power_flow_solver(Ymat,
Pbus, Qbus, E1)

FUNCTION INPUTS:
Ymat - an N*N matrix, defining the admittances of the network.
admittance = 1/impedance). All values are complex (representing
magnitude and angle). Y(i,i) - is the shunt admittance
of bus i to ground. Y(i,j) - is the admittance connecting
bus i to bus j (given i~=j). If no line connects i and j,
the associated admittance is zero.
Mark 1: do not get confused: Y=Inf is a 'short'. Y=0 is an open circuit.
Infinite admittance result in multiple solutions and are forbidden.
Mark 2: Ymat must be symmetric.
Mark 3: this is NOT the Y matrix shown in textbooks,
but a simple definition of the network's admittances. minus signs and summations must not be used. Those are
added automatically in the code. Just enter the simple admittances of the network.

Pbus - Active power vector. Pbus(i) is the active
power injected to bus i. power is positive for generators and
negative for loads. If no source exists, define Pbus(i)=0.
to prevent possible confusion, a power source is not allowed
on bus 1 (the slack bus).

Qbus - Reactive power vector. Same definitions as Pbus.

E1 - The voltage of the 'slack bus' (bus 1). A usual
selection is E1=1.0 . This is a scalar.

FUNCTION OUTPUTS:

Ebus - The resulting vector of bus voltages. Complex
voltage values (magnitude & phase).

Ibus - The resulting currents' vector (complex values).
Ibus(i) is the current supplied to the network by the
source at bus i. For example: Ibus(1) is the current
supplied by the slack bus. The injected power can
be computed by Ebus(i)*conj(Ibus(i)).
If no power source is present at bus i, Ibus(i) should
result in 0. Ibus(i) includes the current supplied to
shunt admittances at bus i. Thus, the shunt admittance
increases this value.

Imat - an N*N matrix of line currents. Imat(i,j)is the
current flowing from bus i to bus j. Imat(i,i) is the
current in the shunt admittance, not including the current
supplied by the source at that node. The sum of the i'th
raw of Imat should equal Ibus(i)

iter - number of iterations done until convergence