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

در آموزش‌های پیشین مجله فرادرس، با مفاهیم و معادلات پخش بار آشنا شدیم. همچنین روش‌های پخش بار نیوتن رافسون و گوس سایدل را به طور مفصل شرح دادیم. پخش بار نیوتن رافسون در متلب نیز یکی دیگر از آموزش‌هایی بود که به بحث درباره آن پرداختیم. در این آموزش، برنامه‌ پخش بار گوس سایدل در متلب را ارائه خواهیم‌ کرد.

قبل از بررسی برنامه‌های متلب، پیشنهاد می‌کنیم آموزش‌های «پخش بار در سیستم قدرت» و «پخش بار گوس سایدل» را مطالعه کنید.

الگوریتم پخش بار گوس سایدل

در این بخش، خلاصه‌ای از مراحل روش گوس سایدل را بیان می‌کنیم. گام‌های روش پخش بار گوس سایدل به صورت زیر هستند:

گام ۱: مقداردهی اولیه $$\bar{V}_j ^{(0)}=V_j^{sp} \angle 0^ \circ$$ را برای $$j=2, 3,  \ldots , m$$، و $$\bar{V}_j ^{(0)}=1.0 \angle 0^ \circ$$ را برای $$j=(m+1), (m+2),  \ldots , n$$ انجام دهید. شمارنده تکرار را $$k=1$$‌ قرار دهید.

فیلم آموزش دینامیک سیستم های قدرت ۱

اینجا کلیک کنید

گام ۲: برای $$i=2, 3, \ldots, m$$، عملیات زیر را انجام دهید:

الف) توان راکتیو را محاسبه کنید:

$$ \large Q _ { i } ^ { ( k ) } = \sum _ { j = 1 } ^ { n } V _ { i } ^ { ( k - 1 ) } V _ { j } ^ { ( k - 1 ) } Y _ { i j } \sin \left ( \theta _ { i } ^ { ( k - 1 ) } - \theta _ { j } ^ { ( k - 1 ) } -\alpha _ { i j } \right ) $$

ب) اگر $$Q_i^{min} \le Q_i ^{(k)} \le Q_i ^{max}$$ برقرار باشد، مقادیر $$\left |\bar{V}_i ^{(k)} \right | = V_i ^{sp}$$ و $$\theta _i ^{(k)}= \angle \left(A_i^{(k)} \right )$$ را در نظر بگیرید.

مقدار $$A_i^{(k)}$$ از رابطه زیر قابل محاسبه است:

$$ \large A _ { i } ^ { ( k ) } = \frac { 1 } { \bar { Y } _ {i i } } \left [ \frac { P _ { i } - j Q _ { i } ^ { ( k ) } } { \left\{ \bar { V } _ { i } ^ { ( k - 1 ) } \right \} ^ { * } } - \sum _ { j = 1 } ^ { i - 1 } \bar { Y } _ { i j } \bar { V } _ { j } ^ { ( k ) } - \sum _ { j = i + 1 } ^ { n } \bar { Y } _ { i j } \bar { V } _ { j } ^ { ( k - 1 ) } \right ] $$

ج) اگر $$Q_i ^{(k)} \ge Q_i^{max}$$، آن‌گاه مقدار ولتاژ زیر را محاسبه کنید:

$$ \large \bar { V } _ { i } ^ { ( k ) } = \frac { 1 } { \bar { Y } _ { i i } } \left [ \frac { P _ { i } - j Q _ { i } ^ { \max } } { \left \{ \bar { V } _ { i } ^ { ( k - 1 ) } \right \} ^ { * } } - \sum _ { j = 1 } ^ { i - 1 } \bar { Y } _ { i j } \bar { V } _{ j } ^ { ( k) } - \sum _ { j = i + 1 } ^ { n } \bar { Y } _ { i j } \bar { V } _ { j } ^ { ( k - 1 ) } \right ] $$

د) اگر $$Q_i ^{(k)} \le Q_i^{min}$$، آن‌گاه ولتاژ را از رابطه زیر به‌دست آورید:

$$ \large \bar { V } _ { i } ^ { ( k ) } = \frac { 1 } { \bar { Y } _ { i i } } \left [ \frac { P _ { i } - j Q _ { i } ^ { \min } }{ \left \{ \bar { V } _ { i } ^ { ( k - 1 ) } \right \} ^ { * } } - \sum _ { j = 1 } ^ { i - 1 } \bar { Y } _ { i j } \bar { V }_ { j } ^ { ( k ) } - \sum _ { j = i + 1 } ^ { n } \bar { Y } _ { i j } \bar { V } _ { j} ^ { ( k - 1 ) } \right ] $$

گام ۳: برای $$i=(m+1), \ldots , n$$، ولتاژ را از رابطه زیر محاسبه کنید:

$$ \large \bar { V } _ { i } ^ { ( k ) } = \frac { 1 } { \bar { Y } _{ i i } } \left [ \frac { P _ { i } - j Q _ { i } ^ { ( k ) } } { \left \{ \bar { V } _ { i } ^ { ( k - 1 ) } \right \} ^ { * } } - \sum _ { j = 1 } ^ { i - 1 } \bar { Y } _ { i j } \bar { V } _ { j } ^ { ( k ) } - \sum _ { j = i + 1 } ^ { n } \bar { Y } _ { i j } \bar { V } _ { j } ^ { ( k - 1 ) } \right ] $$

گام ۴: خطای زیر را برای همه مقادیر $$i=2, \ldots , n$$ محاسبه کنید:

$$ \large e _ { i } ^ { ( k ) } = \left | \bar { V } _ { i } ^ { ( k ) } - \bar { V } _ { i } ^ { ( k - 1 ) } \right | $$

گام 5: از بین خطاهایی که در مرحله قبل محاسبه کرده‌اید، بزرگ‌ترین مقدار را پیدا کنید:

$$ \large e ^ { ( k ) } = \max \left ( e _ { 2 } ^ { ( k ) } , e _ { 3 } ^ { ( k ) } , \cdots \cdots e _ { n } ^ { ( k ) } \right ) $$

پخش بار گوس سایدل در متلب

لازم به ذکر است که این برنامه برای سیستمی با ۴ شین نوشته شده و می‌توانید با تغییر اطلاعات ورودی آن را برای سیستم مورد نظر خود اجرا کنید.

فیلم آموزش بررسی سیستم های قدرت ۱

اینجا کلیک کنید

تابع input_bus_data حاوی اطلاعات مربوط به شین‌ها است که کاربر می‌تواند آن را تغییر دهد. برنامه این تابع به صورت زیر است:

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

end

تابع input_line_data نیز برای ورود اطلاعات مربوط به خطوط است که کاربر آن‌ها را وارد می‌کند. کد متلب زیر، برنامه این تابع را نشان می‌دهد:

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

end

تابع مهم دیگر، Ybus_matrix است که ماتریس ادمیتانس شین را تشکیل می‌دهد. کد متلب این برنامه به صورت زیر است.

function Ybus = Ybus_matrix()
%% From the line data given, generating the Ybus matrix.
linedata = input_line_data();

init_bus = linedata(:,1);       % From bus number
final_bus = linedata(:,2);      % To bus number
r = linedata(:,3);              % Resistance, R
x = linedata(:,4);              % Reactance, X
g = linedata(:,5);              % Conductance, G
b = linedata(:,6);              % Susceptance, B
s = linedata(:,8);              % Shunt or Ground Admittance
z = r + 1i*x;                   % Impedance
y = g + 1i*b;                   % Admittance
s = 1i*s;                       % Shunt from bus to the ground

tot_buses = max(max(init_bus),max(final_bus));    % total no. of buses
tot_branches = length(init_bus);                  % no. of branches
Ybus = zeros(tot_buses,tot_branches);             % Initialising YBus

% Creating a shunt bus to store only the shunt admittances.
sbus = zeros(tot_buses,tot_branches);             % Initialising Shunt bus

% Generating Ybus matrix.
for a = 1:tot_buses
    Ybus(init_bus(a),final_bus(a)) = -y(a);
    Ybus(final_bus(a),init_bus(a)) = -y(a);
end

for b = 1:tot_buses
    sbus(init_bus(b),final_bus(b)) = s(b);
    sbus(final_bus(b),init_bus(b)) = s(b);
end

for b = 1:tot_buses
    for c = 1:tot_buses
        if c ~= b
            Ybus(b,b) = Ybus(b,b) - Ybus(b,c) + sbus(b,c);
        end
    end
end
% Note:- If the Ybus matrix is already known to us then we can comment out
% this portion of the code and directly give Ybus matrix as the input.

end

برنامه اصلی پخش بار گوس سایدل در متلب نیز به صورت زیر است:

% Note:- This code is applicable for any number of buses.
% For an illustration a 4-bus example is considered here.
clear all;
close all;
%% Base Values
V_base = 230e3;     % Base Voltage
MVA_base = 100;     % Base MVA 
alpha = 1.6;        % Accelarating Factor

%% Calling function for obtaining line and bus data.
linedata = input_line_data();
busdata = input_bus_data();
Ybus = Ybus_matrix();

%% Compute reactive power 'Q' for the voltage controlled buses.
bus_no = busdata(:,1);     % Bus Number
bus_type = busdata(:,2);   % Bus Type
GenMW = busdata(:,3);      % Active Power Generated
GenMVAR = busdata(:,4);    % Reactive Power Generated
LoadMW = busdata(:,5);     % Active Power Demanded
LoadMVAR = busdata(:,6);   % Reactive Power Demanded
V = busdata(:,7);          % Initial Bus Voltages
del = busdata(:,8);        % Initial Bus Angles
Qmin = busdata(:,9);       % Minimum limit on the reactive power
Qmax = busdata(:,10);      % Maximum limit on the reactive power

V_orig = V.*cos(del) + 1i*V.*sin(del);
V_new = V_orig;                       % V_orig for storing previous iteration values.
P = (GenMW - LoadMW)/MVA_base;        % Pi = PGi - PLi, Active Power at i'th bus.
Q = (GenMVAR - LoadMVAR)/MVA_base;    % Qi = QGi - QLi, Reactive Power at i'th bus.
        
%% Starting the iterations.
% The iteration will continue until the difference between two consecutive
% voltage values becomes less than 0.00001.

tot_buses = length(bus_no);
iter = 1;
disp("Per unit voltage of buses after each iteration are:")
while(V_orig - V_new < 0.00001)
for idx = 1:tot_buses
    temp1 = 0;
    % Computing new voltages for all the load buses.
    if bus_type(idx) == 3
        for b = 1:tot_buses
            if b ~= idx
                temp1 = temp1 + Ybus(idx,b)*V_new(b);
            end
        end
        Vidx_new = ((P(idx)-1i*Q(idx))/V_new(idx) - temp1)/Ybus(idx,idx);
        Vidx_new_acc = (1-alpha)*V_new(idx) + alpha*Vidx_new;
        V_new(idx) = Vidx_new_acc;
    end
    
    % Computing Q values for all the voltage controlled buses.
    for a = 1:tot_buses
        temp = 0;
        if bus_type(a) == 2
            for b = 1:tot_buses
                temp = temp + V_new(a)*Ybus(a,b)*V_new(b);
            end
            Q(a) = -imag(temp);
            if Q(a) < Qmin(a)
                Q(a) = Qmin(a);
            end
            if Q(a) > Qmax(a)
               Q(a) = Qmax(a);
            end
        end
    end

    % Computing new voltages for the voltage controlled buses.
    if bus_type(idx) == 2
        for b = 1:tot_buses
            if b ~= idx
                temp1 = temp1 + Ybus(idx,b)*V_new(b);
            end
        end
        Vidx_new = ((P(idx)-1i*Q(idx))/V_new(idx) - temp1)/Ybus(idx,idx);
        V_temp = Vidx_new*V_new(idx)/abs(Vidx_new);
        V_new(idx) = V_temp;
    end
end

disp(['Iteration no: ',num2str(iter)])
disp(V_orig)

if max(V_orig - V_new) > 0.00001
    V_orig = V_new;
end
iter = iter + 1;
end

اگر این برنامه را اجرا کنیم، خروجی آن به صورت زیر خواهد بود:

اگر علاقه‌مند به یادگیری مباحث مشابه مطلب بالا هستید، آموزش‌هایی که در ادامه آمده‌اند نیز به شما پیشنهاد می‌شوند:

• مجموعه آموزش‌های مهندسی قدرت
• آموزش بررسی سیستم های قدرت 1
• مجموعه آموزش‌های مهندسی برق
• آموزش نرم‌ افزار DIgSILENT برای آنالیز و شبیه سازی سیستم های قدرت
• پایداری سیستم قدرت — به زبان ساده
• شبکه عصبی در متلب — از صفر تا صد
• پارامترهای خط انتقال در مهندسی قدرت — به زبان ساده

^^