MATLAB Coding for DC Power Flow

Let us consider that the conductance ‘G’ of a transmission line is very low as compared to the susceptance ‘B’. Also under normal operating conditions of the transmission line, the difference in the angles of the voltage at the two buses, which are connected by the given transmission line, is usually very less (less than 15^o). Similarly, the numerical values of the voltage at the two buses are very close to 1.0 p.u. With these approximations, the real power flow in a transmission line is proportional to the circuit susceptance ‘B’ and the difference in voltage phasor angles.

The DC power flow equations (in matrix form) based on the above fact is given as:

[P] = [B’][theta]

where [P] is the vector of nodal active power injection for all the buses except the reference bus,

[theta] is the vector of nodal phase angles for all the buses except the reference bus,

[B’] is the B-prime matrix

Now [B’]= [D] *[ S]

where [D] is a matrix having non-diagonal elements of zeros and the diagonal elements in position row m, column m contains the negative of the susceptance of the m^th branch.

S is branch-node incidence matrix also called the adjacency matrix or the connection matrix after eliminating a column corresponding to reference bus.

The MATLAB Coding for the DC power flow is as given below:

The inputs required for the coding are the branch data and the bus data. Branch data should have the numbering of the total branches in the network, from bus number, to bus number, and branch susceptance. The bus data should have the numbering of each buses in the network, and load and generation data at each bus.

% “br_data” is the branch data; “bus” is the bus data.

% Col.1 of “br_data” is the numbering of branches.

% Col.2 of “br_data” is the ‘from-bus’ number.

% Col.3 of “br_data” is the ‘to-bus’ number.

% Col.4 of “br_data” is the susceptance of branches.

% Col.1 of “bus” is the numbering of buses.

% Col.2 of “bus” is the active load at respective buses.

% Col.3 of “bus” is the active generation at respective buses.

% mention the reference bus number

>> ref_bus=6; % here bus number 6 is considered the reference bus

% Y-BUS formation

>>  YBUS=Y_bus(br_data);

% function file ‘Y_bus’is called upon; please refer to previous blogs on related topic

>> B= imag(YBUS);

% Eliminating the corresponding row & col. i.e. row & col. of reference bus

>> B(ref_bus, :)=[];

>> B(:, ref_bus,)=[];

>> BB= -B;

% Calculate the active power injected, i.e. the difference in generation and the demand at each bus

>> Pbus=(bus(:, 3) - bus(:, 2));

% Calculate the voltage angle

>> theta=(BB)^ -1* Pbus;

>> y = imag(br_data(:, 4))’;

>> D = -diag(y);

% formation of “bus-incidence” matrix

>> elements = max(br_data(:, 1));

>> for i =1:elements,

>> if br_data(i,2) ~= 0,

>> S(i, br_data(i,2)) = 1;

>> end

>> if br_data(i,3) ~= 0,

>> S(i, br_data(i,3)) = -1;

>> end

>> end

% Calculation for line flows (active power flows)

>> P = D*S*theta