Gauss Seidel method for power flow analysis

Question

I need to write a matlab code for a power flow analysis using the Gauss-Seidel method. I have done the first part of the question, but I cannot figure out how to do the second one. The figure given on the assignment is the following:

The question asks for a code to be written for a power mismatch tolerance of Ep=10^(-6) and max(|ΔP2|,|ΔQ2|)<Ep.

The code I've written so far is the following:

% In the following code the Gauss-Seidel method will be used in both parts
% (1) and (2) to calculate the power flow in a two-bus system.
% Bus 1 is a slack bus, for which the voltage and angle are given.
% Bus 2 is a PQ bus, for which the load power is given, but the voltage and
% angle need to be computed.


% Part (1):
% In this section, the power flow solution to be found has a voltage
% convergence tolerance of Ep =1e-6, with max |Vi(k+1)-Vi(k)|< Ep.

% Formation of the matrix for the bus admittance:
% The transmition impedance between buses 1 and 2 is z12 = 0.0 + j0.12 p.u. 
% therefore: 
z12=0.0+1i*0.12;
% y12 = 1 / z12
y12=1/z12;
Y11=y12;
Y12=-y12;
Y21=Y12;
Y22=y12;

% Formation of power injections:
PQ=[]; % the bus number is represented by the number in the brackets, i.e. 
       % PQ(1) is for bus 1, and PQ(2) is for bus 2
% The load at bus 2 is 0.90 + j0.25 p.u., thus the power injection is:
PQ(2)=-0.90-1i*0.25;

% The follwoing are the symbols of the power flow solutions' values:
VK=[];
VK1=[];
IK=[];

% The voltage at bus 1 is:
VK(1)=1.05+1i*0.0;
VK1(1)=1.05+1i*0.0;

% The voltage at bus 2, where the flat start is, is the following:
VK(2)=1.0+1i*0.0;
VK1(2)=1-0+1i*0.0;
% Also, the current at bus 2 starts at:
IK(2)=0.0+1i*0.0;

%Convergence tolerance format:
EP=input('Convergence Tolerance, EP:');
fprintf('\n');
fprintf('%12.5e\n', EP);

%Following are the Gauss-Seidel iterations:
for k=1:1000
% The current at bus 2 is calculated according to the equation: IK(2) =
% conjugate ((P(2) + jQ(2)/VK(2)), so:
IK(2)=conj(PQ(2)/VK(2));
% The voltage at bus 2 (which is the PQ bus), is given by the equation:
% VK1(2) = VK+1(2)/Y22, where VK+1(2) = (IK(2)-Y21xV1,therefore:
VK1(2)=(IK(2)-Y21*VK1(1))/Y22; 
fprintf('The number of iterations is the following: %4i max|VK+1-VK|=%12.5e.\n The voltage magnitude at bus 2 is: %12.9e.\n The voltage angle at bus 2 is equal to: %12.9e.\n',k,abs(VK1(2)-VK(2)),abs(VK1(2)),angle(VK1(2))*180/pi);
% Convergence calculation:
    if (abs(VK1(2)-VK(2))>=EP)
    % In the case max|VK+1-VK|>=EP, then stop:
        break;
    else
    % In any other situation, the iteration should be continued, and the
    % solution, for which VK=VK+1, should be updated:
        VK(2)=VK1(2);
    end
end

% The current injection at bus 1 (slack bus) is:
IK(1)=Y11*VK1(1)+Y12*VK1(2);
% The power injection at bus 1 is determined using the form:
PQ(1)=VK1(1)*conj(IK(1));

% Printing the output power injections (active and reactive), at slack
% bus 2:
fprintf('Reactive power injection at bus 2: %12.5e. Active power injection at bus 2: %12.5e\n', real(PQ(1)), imag(PQ(1)));


% Part(2):
% In this section, the power flow solution to be found has a voltage
% convergence tolerance of Ep =1e-6, with max (|ΔP2|,|ΔQ2|)< Ep.

I am new to power flow analysis and matlab so helping me with this would be extremely useful. Is there any part of the code that I've written so far that I should use for the second part, and also, how will I calculate ΔP2 and ΔQ2. Thanks in advance.

Answer 1

My expectation is that the second part is just requiring you to go beyond only using change in bus voltage magnitude at each step (to determine when the case is solved) and adding a check for the \$P\$ and \$Q\$ tolerance. I would think you would want to calculate the total \$P\$ and \$Q\$ at your bus 2 and sum it with the load injection at that bus. If you have a perfect solution the result will be zero. If result > \$\epsilon\$, iterate again.

Using this figure as a reference,

The power injection at your bus i will be,

$$P=\sum_{j=1}^n|V_i||Y_{ij}||V_j|cos(\delta _j+\theta_{ij}-\delta_i) $$

From this you subtract the scheduled load for that bus, \$P_i\$, and the closer to zero the better the solution. This is what I think your \$\Delta P\$ will be.

Do similar for reactive with,

$$Q=\sum_{j=1}^n|V_i||Y_{ij}||V_j|sin(\delta _j+\theta_{ij}-\delta_i)$$

I highly recommend you use another tool to check yourself. MATPOWER is free and runs on top of Matlab (or Octave if you prefer a free equivalent to Matlab). I created a test case file like yours (I named it test.m) and set the bus voltage base to 1.0kV to leave it in per-unit.

function mpc = test
%% MATPOWER Case Format : Version 2
mpc.version = '2';
%%-----  Power Flow Data  -----%%
%% system MVA base
% This is the base that the branch impedances are on.  
mpc.baseMVA = 100;
%% Vm under mpc.gen sets the slack bus (1) voltage magnitude
%% Va under mpc.bus sets the slack bus (1) voltage angle 
%   bus_i   type      Pd       Qd      Gs       Bs  area      Vm      Va      baseKV    zone     Vmax      Vmin
mpc.bus = [
  1    3     0      0     0      0     1     1     0      1.0     1      1       1;
  2    1  0.90   0.25     0      0     0     1     1      1.0     1      1       0;
];
%% generator data
%   bus    Pg       Qg    Qmax    Qmin        Vg        mBase   status  Pmax    Pmin    Pc1 Pc2 Qc1min  Qc1max  Qc2min  Qc2max  ramp_agc    ramp_10 ramp_30 ramp_q  apf   ?
mpc.gen = [ 
     1      0       0        0       0      1.05     0       1    0     0    0   0       0       0       0       0         0        0       0      0    0   0;
];
%% branch data
%   fbus    tbus             r           x      b      rateA    rateB   rateC    ratio   angle   status   angmin      angmax
mpc.branch = [
 1     2         0.0        0.12   0        0      0     0      1        0        1        0         0;
];

I ran it with the runpf ('test') command in Octave and the results from Newton Rahpson load flow was bus 2 voltage of 1.05∠-0.056° pu. I didn't try real hard to check the data, I just wanted to point you to a way to check your own power flow results.

Also, be careful with your code notes, could trip you up later. For example you have:

% conjugate ((P(2) + jQ(2)/VK(2)), so:

But I think you mean, "% conjugate [(P(2) + jQ(2)]/VK(2), so:"