پخش بار گوس سایدل در متلب — راهنمای کاربردی
در آموزشهای پیشین مجله فرادرس، با مفاهیم و معادلات پخش بار آشنا شدیم. همچنین روشهای پخش بار نیوتن رافسون و گوس سایدل را به طور مفصل شرح دادیم. پخش بار نیوتن رافسون در متلب نیز یکی دیگر از آموزشهایی بود که به بحث درباره آن پرداختیم. در این آموزش، برنامه پخش بار گوس سایدل در متلب را ارائه خواهیم کرد.
قبل از بررسی برنامههای متلب، پیشنهاد میکنیم آموزشهای «پخش بار در سیستم قدرت» و «پخش بار گوس سایدل» را مطالعه کنید.
الگوریتم پخش بار گوس سایدل
در این بخش، خلاصهای از مراحل روش گوس سایدل را بیان میکنیم. گامهای روش پخش بار گوس سایدل به صورت زیر هستند:
گام ۱: مقداردهی اولیه $$\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 حاوی اطلاعات مربوط به شینها است که کاربر میتواند آن را تغییر دهد. برنامه این تابع به صورت زیر است:
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
تابع input_line_data نیز برای ورود اطلاعات مربوط به خطوط است که کاربر آنها را وارد میکند. کد متلب زیر، برنامه این تابع را نشان میدهد:
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
تابع مهم دیگر، Ybus_matrix است که ماتریس ادمیتانس شین را تشکیل میدهد. کد متلب این برنامه به صورت زیر است.
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
اگر این برنامه را اجرا کنیم، خروجی آن به صورت زیر خواهد بود:
اگر علاقهمند به یادگیری مباحث مشابه مطلب بالا هستید، آموزشهایی که در ادامه آمدهاند نیز به شما پیشنهاد میشوند:
- مجموعه آموزشهای مهندسی قدرت
- آموزش بررسی سیستم های قدرت 1
- مجموعه آموزشهای مهندسی برق
- آموزش نرم افزار DIgSILENT برای آنالیز و شبیه سازی سیستم های قدرت
- پایداری سیستم قدرت — به زبان ساده
- شبکه عصبی در متلب — از صفر تا صد
- پارامترهای خط انتقال در مهندسی قدرت — به زبان ساده
^^
سلام ممنون از مطلب مفیدتون،
یه سوال
در تابع input_Line_data
اون مقدار Total Charging
چطور به دست میاد ؟
Pgen-Pd هست آیا ؟ بر حسب پریونیته یا چیز دیگه ایه ؟