Politechnika Częstochowska
Wydział Elektryczny
Magisterskie Studia Uzupełniające
LABORATORIUM
GOSPODARKA ELEKTROENERGETYCZNA
Ćwiczenie: Zastosowanie procedur programowania liniowego oraz optymalizacji nieliniowej w pakiecie MatLab do rozwiązywania zadań z optymalizacji gospodarki elektroenergetycznej
Wykonał:
Krystek Szymon
1. Cel ćwiczenia
Celem ćwiczenia jest zapoznanie się z podstawowymi zasadami programowania i rozwiązywania zadań z elektroenergetyki za pomocą programu Matlab. Aby rozwiązać zadanie w Matlabie należy stworzyć odpowiednie funkcje i skrypty.
2. Zasady programowania
2.1 Tworzenie skryptów
Skrypt jest plikiem tekstowym zawierającym zestaw funkcji i poleceń Matlaba. Pliki skryptowe mają rozszerzenie .m. Pliki skryptowe można tworzyć w każdym edytorze tekstowym. Najwygodniej wykorzystać edytor Matlaba. Dostęp do edytora jest możliwy przez odpowiednią ikonę.
2.2 Opisywanie skryptów
Każdy skrypt powinien mieć krótki opis zawartości i działania. Opis umieszcza się za
znakiem %.
2.3 Tworzenie funkcji
Funkcja tak jak skrypt jest plikiem tekstowym zawierającym zestaw funkcji i poleceń, Matlaba i zaczynać się powinna od słowa kluczowego function. Pliki funkcji mają również rozszerzenie .m. Ważne jest aby nazwa funkcji i nazwa pliku były takie same. Pliki funkcji można tworzyć w każdym edytorze tekstowym. Najwygodniej wykorzystać edytor Matlaba. Podstawową różnicą miedzy funkcją a skryptem jest sposób przechowywania danych. Skrypt
czyni to w przestrzeni roboczej, natomiast funkcja przechowuje je poza przestrzenia roboczą,
co pozwala na dublowanie nazw zmiennych z przestrzenią roboczą. Inaczej mówiąc funkcja jest hermetyczna i pokazuje na zewnątrz tylko dane wyjściowe.
Na przykładzie takiego programowania zostało rozwiązane zadanie 1 i zadanie 2.
Zadanie 1
Trzy elektrownie należące do jednego koncernu energetycznego zasilają węzeł sieciowy (GPZ Główny Punkt Zasilający). Chwilowe zapotrzebowanie na moc w GPZ wynosi (200+4*15)=260 MW przy cosϕ=0,8. Zakładając, że każda elektrownia może pokryć powyższe zapotrzebowanie, wyznaczyć wartości mocy dosłanych do węzła GPZ z elektrowni E1, E2 i E3, minimalizujących straty mocy czynnej w sieci przesyłowej. Dane linii przesyłowych:
L1: E1-GPZ Un=110 kV, l1=84 km, S1=240 mm2, γ=35 MS/m,
L2: E2-GPZ Un=110 kV, l2=84 km, S2=120 mm2, γ=35 MS/m,
L3: E3-GPZ Un=110 kV, l3=136,5 km, S3=520 mm2, γ=35 MS/m,
Do rozwiązania tego zadania trzeba napisać skrypty:
function dP=mojafun(S)
dP=( 1/110)^2*(10*S(1)^2+20*S(2)^2+7.5*S(3)^2);
Następnie
function [c,ceq]=mojeogr(S)
c=[];
ceq=[S(1)+S(2)+S(3)-236/0.8];
Oraz na końcu skrypt wywołujący
%skrypt wywołujący
S0 = [50;50;50];
lb = [0;0;0];
ub = [120;80;200];
[S,fval,exitflag]= fmincon('mojafun',S0,[],[],[],[],lb,ub,'mojeogr')
%koniec
Rozwiązanie zadania:
Active Constraints:
1
S =
104.1177
52.0583
138.8240
fval =
25.3841
exitflag =
1
Zadanie 2
W stacji transformatorowej 15/0,4 kV/kV zainstalowane są trzy transformatory:
T1: SNT1=400 kVA, ΔPj=1,3 kW, ΔPobc N=4,800 kW
T2: SNT2=630 kVA, ΔPj=1,8 kW, ΔPobc N=7,938 kW
T3: SNT3=1000 kVA, ΔPj=2,6 kW, ΔPobc N=15,000 kW
Dla chwilowego obciążenia w stacji po stronie 0,4 kV Pobc= (1,62+0,03*15)=2,07 MW przy cosϕ=0,9 wyznaczyć obciążenie każdego z transformatorów przy założeniu minimum strat mocy czynnej w transformatorach.
Do rozwiązania tego zadania trzeba napisać skrypty:
Pierwszym skryptem będzie skrypt przedstawiony kodem
function deltaP=stratywtrafo(S)
global w1 w2 w3
deltaP=w1*(1.3+4.8*(S(1)/0.4).^2)+w2*(1.8+7.938*(S(2)/0.63).^2)+w3*(2.6+15*(S(3)/1).^2);
Następnie trzeba napisać skrypt wywołujący
%skrypt wywpłujący
global w1 w2 w3
w1=0
w2=0
w3=0
Pobc=1.89
cosfi=0.9
Sobc=Pobc/cosfi
if Sobc>1.16897
w1=1;
w2=1;
w3=1;
elseif Sobc>0.94083
w2=1;
w3=1;
elseif Sobc>0.64795
w1=1;
w3=1;
elseif Sobc>0.40875
w1=1;
w2=1;
elseif Sobc>0.400
w3=1;
elseif Sobc>0.22360
w2=1;
else
w1=1;
end
Aeq=[w1 w2 w3];
beq=[Sobc];
lb=[0;0;0;];
ub=[0.48;0.756;1.2]
S0=[0;0;0];
[S,deltaP,exitflag]=fmincon('stratywtrafo',S0,[],[],Aeq,beq,lb,ub)
W zadaniu 2 „w” oznaczają pracę odpowiednich transformatorów przy danym obciążeniu. Obciążenie chwilowe w stacji po stronie 0,4kV w moim przypadku wyniosło 2,07 MW.
Otrzymanym wynikiem będzie:
w1 =
0
w2 =
0
w3 =
0
Pobc =
1.8900
cosfi =
0.9000
Sobc =
2.1000
ub =
0.4800
0.7560
1.2000
Optimization terminated successfully:
Magnitude of directional derivative in search direction
less than 2*options.TolFun and maximum constraint violation
is less than options.TolCon
Active Constraints:
1
5
6
S =
0.4800
0.7560
0.8640
deltaP =
35.2402
exitflag =
1
Zadanie 3
Ze stacji RG 110/6 kV/kV zasilane są trzy rozdzielnie oddziałowe RO1, RO2 i RO3 następującymi kablami 6 kV:
Tor RG-RO1: YAKYy 3x240, l1=800 m
Tor RG-RO2: YAKYy 3x120, l2=200 m
Tor RG-RO3: YAKYy 3x50, l3=500 m
Dla chwilowego obciążenia poszczególnych rozdzielni:
RO1: P1=1 600 kW, cos ϕ1 = 0,8
RO2: P2=1 000 kW, cos ϕ2 = 0,85
RO3: P3=630 kW, cos ϕ3 = 0,8
Dobrać baterię kondensatorów potrzebnych do utrzymania w RG współczynnika mocy tg ϕk=0,4. Rozdzielić baterię kondensatorów do kompensacji grupowej, przy założeniu minimum strat mocy czynnej na terenie przedsiębiorstwa.
Rozwiązaniem podobnie do wcześniejszych zadań piszemy skrypty:
function deltaPq=kompgrupowa(Qki)
global U Qi Ri
deltaPq=(1/(U^2))*(((Qi(1)Qki(1))^2*Ri(1))+((Qki(2))^2*Ri(2))+((Qi(3))^2*Ri(3)));
następnie piszemy skrypt wywołujący
%skrypt wywołujący kompensacja grupowa
global U Ri Qi
U=6
Pi=[1600 1600 630]
cos=[0.8 0.85 0.8]
Qi=((Pi./cos).^2-Pi.^2).^(1/2)
disp ('Moce bierne')
Si=[240 240 50]
li=[1000 300 600]
g=35;
disp ('Rezystancje linii')
Ri=1/g*li./Si
[m n]=size(Pi);
P=0;
Q=0;
for i=1:n
P=P+Pi(i);
Q=Q+Qi(i);
end
disp ('Tangens przed kompensacją')
tgfi=Q/P
tgfik=0.4
disp ('Moc baterii kondensatorów w kVar')
Qk=P*(tgfi-tgfik)
Q0=[0 0 0]
Aeq=[1 1 1]
beq=[Qk]
lb=[0 0 0]
disp ('Moc baterii kondensatorów w rozdzielniach oddziałowych RO w kVar')
[Qki, deltaPq, exitflag]=fmincon('kompgrupowa',Q0,[],[],Aeq,beq,lb,[])
Wynikiem będzie:
U =
6
Pi =
1600 1600 630
cos =
0.8000 0.8500 0.8000
Qi =
1.0e+003 *
1.2000 0.9916 0.4725
Moce bierne
Si =
240 240 50
li =
1000 300 600
Rezystancje linii
Ri =
0.1190 0.0357 0.3429
Tangens przed kompensacją
tgfi =
0.6956
tgfik =
0.4000
Moc baterii kondensatorów w kVar
Qk =
1.1321e+003
Q0 =
0 0 0
Aeq =
1 1 1
beq =
1.1321e+003
lb =
0 0 0
Moc baterii kondensatorów w rozdzielniach oddziałowych RO w kVar
Active Constraints:
1
3
4
Qki =
1.0e+003 *
1.1321 0 0.0000
deltaPq =
2.1415e+003
exitflag =
1
9