پخش بار گوس سایدل در متلب — راهنمای کاربردی

1function busdata = input_bus_data()
2%% Bus data to be given by the user as the input.
3%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4% Note:-    Bus Types       -      Symbol
5%           Slack Bus       -        1
6%         Generator Bus     -        2
7%           Load Bus        -        3
8%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
9%                   |  Generation |      Load     |
10%          |Bus|Type|P(MW)|Q(MVAR)|P(MW)| Q(MVAR) | V,pu |delta|Qmin|Qmax|
11busdata = [  1   1     0      0      50    30.99    1.00   0.0    0   0.0 ;
12             2   3     0      0     170   105.35    1.00   0.0    0   0.0 ;
13             3   3     0      0     200   123.94    1.00   0.0    0   0.0 ;
14             4   2    318     0      80    49.58    1.02   0.0   0.1  0.2;];
15
16end

1function linedata = input_line_data()
2%% Line data to be given by the user as the input.
3%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4% R - Resistance
5% X - Reactance
6% G - Conductance
7% B - Susceptance
8% Y/2 - Shunt admittance
9%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10%          |From|To |   R    |   X    |   G    |   B    |  Total  |  Y/2  |
11%          |Bus |Bus|   pu   |   pu   |   pu   |   pu   | charging|  pu   |
12linedata = [ 1    2   0.01008  0.05040  3.81563 -19.0781   10.25   0.05125;
13             1    3   0.00744  0.03720  5.16956 -25.8478    7.75   0.03875;
14             2    4   0.00744  0.03720  5.16956 -25.8478    7.75   0.03875;
15             3    4   0.01272  0.06360  3.02371 -15.1185   12.75   0.06375;];
16
17end

1function Ybus = Ybus_matrix()
2%% From the line data given, generating the Ybus matrix.
3linedata = input_line_data();
4
5init_bus = linedata(:,1);       % From bus number
6final_bus = linedata(:,2);      % To bus number
7r = linedata(:,3);              % Resistance, R
8x = linedata(:,4);              % Reactance, X
9g = linedata(:,5);              % Conductance, G
10b = linedata(:,6);              % Susceptance, B
11s = linedata(:,8);              % Shunt or Ground Admittance
12z = r + 1i*x;                   % Impedance
13y = g + 1i*b;                   % Admittance
14s = 1i*s;                       % Shunt from bus to the ground
15
16tot_buses = max(max(init_bus),max(final_bus));    % total no. of buses
17tot_branches = length(init_bus);                  % no. of branches
18Ybus = zeros(tot_buses,tot_branches);             % Initialising YBus
19
20% Creating a shunt bus to store only the shunt admittances.
21sbus = zeros(tot_buses,tot_branches);             % Initialising Shunt bus
22
23% Generating Ybus matrix.
24for a = 1:tot_buses
25    Ybus(init_bus(a),final_bus(a)) = -y(a);
26    Ybus(final_bus(a),init_bus(a)) = -y(a);
27end
28
29for b = 1:tot_buses
30    sbus(init_bus(b),final_bus(b)) = s(b);
31    sbus(final_bus(b),init_bus(b)) = s(b);
32end
33
34for b = 1:tot_buses
35    for c = 1:tot_buses
36        if c ~= b
37            Ybus(b,b) = Ybus(b,b) - Ybus(b,c) + sbus(b,c);
38        end
39    end
40end
41% Note:- If the Ybus matrix is already known to us then we can comment out
42% this portion of the code and directly give Ybus matrix as the input.
43
44end

1% Note:- This code is applicable for any number of buses.
2% For an illustration a 4-bus example is considered here.
3clear all;
4close all;
5%% Base Values
6V_base = 230e3;     % Base Voltage
7MVA_base = 100;     % Base MVA 
8alpha = 1.6;        % Accelarating Factor
9
10%% Calling function for obtaining line and bus data.
11linedata = input_line_data();
12busdata = input_bus_data();
13Ybus = Ybus_matrix();
14
15%% Compute reactive power 'Q' for the voltage controlled buses.
16bus_no = busdata(:,1);     % Bus Number
17bus_type = busdata(:,2);   % Bus Type
18GenMW = busdata(:,3);      % Active Power Generated
19GenMVAR = busdata(:,4);    % Reactive Power Generated
20LoadMW = busdata(:,5);     % Active Power Demanded
21LoadMVAR = busdata(:,6);   % Reactive Power Demanded
22V = busdata(:,7);          % Initial Bus Voltages
23del = busdata(:,8);        % Initial Bus Angles
24Qmin = busdata(:,9);       % Minimum limit on the reactive power
25Qmax = busdata(:,10);      % Maximum limit on the reactive power
26
27V_orig = V.*cos(del) + 1i*V.*sin(del);
28V_new = V_orig;                       % V_orig for storing previous iteration values.
29P = (GenMW - LoadMW)/MVA_base;        % Pi = PGi - PLi, Active Power at i'th bus.
30Q = (GenMVAR - LoadMVAR)/MVA_base;    % Qi = QGi - QLi, Reactive Power at i'th bus.
31        
32%% Starting the iterations.
33% The iteration will continue until the difference between two consecutive
34% voltage values becomes less than 0.00001.
35
36tot_buses = length(bus_no);
37iter = 1;
38disp("Per unit voltage of buses after each iteration are:")
39while(V_orig - V_new < 0.00001)
40for idx = 1:tot_buses
41    temp1 = 0;
42    % Computing new voltages for all the load buses.
43    if bus_type(idx) == 3
44        for b = 1:tot_buses
45            if b ~= idx
46                temp1 = temp1 + Ybus(idx,b)*V_new(b);
47            end
48        end
49        Vidx_new = ((P(idx)-1i*Q(idx))/V_new(idx) - temp1)/Ybus(idx,idx);
50        Vidx_new_acc = (1-alpha)*V_new(idx) + alpha*Vidx_new;
51        V_new(idx) = Vidx_new_acc;
52    end
53    
54    % Computing Q values for all the voltage controlled buses.
55    for a = 1:tot_buses
56        temp = 0;
57        if bus_type(a) == 2
58            for b = 1:tot_buses
59                temp = temp + V_new(a)*Ybus(a,b)*V_new(b);
60            end
61            Q(a) = -imag(temp);
62            if Q(a) < Qmin(a)
63                Q(a) = Qmin(a);
64            end
65            if Q(a) > Qmax(a)
66               Q(a) = Qmax(a);
67            end
68        end
69    end
70
71    % Computing new voltages for the voltage controlled buses.
72    if bus_type(idx) == 2
73        for b = 1:tot_buses
74            if b ~= idx
75                temp1 = temp1 + Ybus(idx,b)*V_new(b);
76            end
77        end
78        Vidx_new = ((P(idx)-1i*Q(idx))/V_new(idx) - temp1)/Ybus(idx,idx);
79        V_temp = Vidx_new*V_new(idx)/abs(Vidx_new);
80        V_new(idx) = V_temp;
81    end
82end
83
84disp(['Iteration no: ',num2str(iter)])
85disp(V_orig)
86
87if max(V_orig - V_new) > 0.00001
88    V_orig = V_new;
89end
90iter = iter + 1;
91end