第一篇：電力系統(tǒng)潮流計算的MATLAB輔助程序設(shè)計,潮流計算程序
電力系統(tǒng)潮流計算的MATLAB輔助程序設(shè)計
潮流計算，通常指負(fù)荷潮流，是電力系統(tǒng)分析和設(shè)計的主要組成部分，對系統(tǒng)規(guī)劃、安全運行、經(jīng)濟調(diào)度和電力公司的功率交換非常重要。此外，潮流計算還是其它電力系統(tǒng)分析的基礎(chǔ)，比如暫態(tài)穩(wěn)定，突發(fā)事件處理等。現(xiàn)代電力系統(tǒng)潮流計算的方法主要:高斯法、牛頓法、快速解耦法和MATLAB的M語言編寫的MATPOWER4.1，這里主要介紹高斯法、牛頓法和快速解耦法。高斯法的程序是lfgauss，其與lfybus、busout和lineflow程序聯(lián)合使用求解潮流功率。lfybus、busout和lineflow程序也可與牛頓法的lfnewton程序和快速解耦法的decouple程序聯(lián)合使用。(讀者可以到MATPOWER主頁下載MATPOWER4.1,然后將其解壓到MATLAB目錄下，即可使用該軟件進行潮流計算)
一、高斯-賽德爾法潮流計算使用的程序：
高斯-賽德法的具體使用方法讀者可參考后面的實例，這里僅介紹各程序的編寫格式: lfgauss：該程序是用高斯法對實際電力系統(tǒng)進行潮流計算，需要用到busdata和linedata兩個文件。程序設(shè)計為輸入負(fù)荷和發(fā)電機的有功MW和無功Mvar，以及節(jié)點電壓標(biāo)幺值和相角的角度值。根據(jù)所選復(fù)功率為基準(zhǔn)值將負(fù)荷和發(fā)電機的功率轉(zhuǎn)換為標(biāo)幺值。對于PV節(jié)點，如發(fā)電機節(jié)點，要提供一個無功功率限定值。當(dāng)給定電壓過高或過低時，無功功率可能超出功率限定值。在幾次迭代之后(高斯-塞德爾迭代為10次)，需要檢查一次發(fā)電機節(jié)點的無功出力，如果接近限定值，電壓幅值進行上下5%的調(diào)整，使得無功保持在限定值內(nèi)。
lfybus：這個程序需要輸入線路參數(shù)、變壓器參數(shù)以及變壓器分接頭參數(shù)。并將這些參數(shù)放在名為linedata的文件中。這個程序?qū)⒆杩罐D(zhuǎn)換為導(dǎo)納，并得到節(jié)點導(dǎo)納矩陣。
busout：該程序以表格形式輸出結(jié)果，節(jié)點輸出包括電壓幅值和相角，發(fā)電機和負(fù)荷的有功和無功功率，以及并聯(lián)電容器或電抗器的有功和無功功率。
lineflow：該程序輸出線路的相關(guān)數(shù)據(jù)，程序設(shè)計輸出流入線路終端的有功和無功的功率、線損以及節(jié)點功率，還包含整個系統(tǒng)的有功和無功損耗。
lfnewton是牛頓-拉夫遜法對實際電力系統(tǒng)潮流計算開發(fā)的程序，數(shù)據(jù)準(zhǔn)備和程序格式和高斯-賽德爾法一樣，包括程序lfybus，busout和lineflow。
decouple是快速解耦法對實際電力系統(tǒng)潮流計算開發(fā)的程序，同高斯法和牛頓法一樣需要用到三個程序：lfybus、busout、lineflow。
二、數(shù)據(jù)準(zhǔn)備
為了在MATLAB環(huán)境下用高斯法進行潮流計算，必須定義下列變量：基準(zhǔn)功率，功率允許誤差，加速因子和最大迭代次數(shù)。上述變量命名(小寫字母)為：basemva、accuracy、accel和maxiter，一般規(guī)定為：basemva=100； accuracy=0.001；accel=1.6；maxiter=80；輸入文件準(zhǔn)備的第一步是給節(jié)點編號，節(jié)點號碼必須是連續(xù)的，但節(jié)點數(shù)據(jù)輸入不一定按順序來編寫。此外，還需要下列數(shù)據(jù)文件：
1.節(jié)點數(shù)據(jù)文件busdata：節(jié)點信息輸入格式為單行輸入，輸入的數(shù)據(jù)形成一個矩陣，叫做busdata矩陣。第一列為節(jié)點號；第二列為節(jié)點類型；第三列和第四列分別為節(jié)點電壓幅值(標(biāo)幺值)和相角(單位為度)；第五列和第六列分別為負(fù)荷的有功功率和無功功率；第七列到十列分別為發(fā)電機的有功功率、無功功率、最小無功出力和最大無功出力；最后一列為并聯(lián)電容器注入無功功率。第二列的編碼用0、1、2來區(qū)分PQ節(jié)點、平衡節(jié)點和PV節(jié)點：
0表示PQ節(jié)點，輸入正的有功功率（MW）和無功功率（Mvar），并且要設(shè)定節(jié)點電壓初始估計值，一般幅值和相角分別設(shè)為1和0，若已經(jīng)給定初始值，則用其給定值來代替1和0。
1表示平衡節(jié)點，且已知該節(jié)點的電壓幅值和相角。
2表示PV節(jié)點，要設(shè)定該節(jié)點的節(jié)點電壓幅值和發(fā)電機的有功功率（MW），并設(shè)定發(fā)電機的無功最小出力和最大出力（Mvar）。
2.線路數(shù)據(jù)文件linedata線路數(shù)據(jù)用節(jié)點對的方法來確定，數(shù)據(jù)包含在稱為linedata的矩陣中。第一列和第二列為節(jié)點號碼，第三列到第五列為線路電阻、電抗及該線路電納值的一半，以標(biāo)幺值表示。最后一列為變壓器分接頭設(shè)定值，對線路來說，需要輸入1。線路輸入為無輸入順序，對變壓器來說，左側(cè)的節(jié)點號設(shè)為分接頭端。
3.zdata是線路數(shù)據(jù)輸入變量，包括四項，前兩項是節(jié)點編號，后兩項是線路電阻和電抗，均以標(biāo)幺值表示，函數(shù)返回節(jié)點導(dǎo)納矩陣。
三、潮流計算的MATLAB程序清單
1.lfgauss.m程序清單
% Power flow solution by Gauss-Seidel method Vm=0;delta=0;yload=0;deltad =0;nbus = length(busdata(:,1));kb=[];Vm=[];delta=[];Pd=[];Qd=[];Pg=[];Qg=[];Qmin=[];Qmax=[];Pk=[];P=[];Qk=[];Q=[];S=[];V=[];for k=1:nbus n=busdata(k,1);kb(n)=busdata(k,2);Vm(n)=busdata(k,3);delta(n)=busdata(k, 4);Pd(n)=busdata(k,5);Qd(n)=busdata(k,6);Pg(n)=busdata(k,7);Qg(n)= busdata(k,8);Qmin(n)=busdata(k, 9);Qmax(n)=busdata(k, 10);Qsh(n)=busdata(k, 11);if Vm(n)<= 0 Vm(n)= 1.0;V(n)= 1 + j*0;else delta(n)= pi/180*delta(n);V(n)= Vm(n)*(cos(delta(n))+ j*sin(delta(n)));P(n)=(Pg(n)-Pd(n))/basemva;Q(n)=(Qg(n)-Qd(n)+ Qsh(n))/basemva;S(n)= P(n)+ j*Q(n);end
DV(n)=0;end
num = 0;AcurBus = 0;converge = 1;Vc = zeros(nbus,1)+j*zeros(nbus,1);Sc = zeros(nbus,1)+j*zeros(nbus,1);
while exist('accel')~=1 accel = 1.3;end
while exist('accuracy')~=1 accuracy = 0.001;end
while exist('basemva')~=1 basemva= 100;end
while exist('maxiter')~=1 maxiter = 100;end
mline=ones(nbr,1);for k=1:nbr for m=k+1:nbr if((nl(k)==nl(m))&(nr(k)==nr(m)));mline(m)=2;elseif((nl(k)==nr(m))&(nr(k)==nl(m)));mline(m)=2;else, end end end
iter=0;maxerror=10;while maxerror >= accuracy & iter <= maxiter iter=iter+1;for n = 1:nbus;YV = 0+j*0;for L = 1:nbr;if(nl(L)== n & mline(L)== 1), k=nr(L);YV = YV + Ybus(n,k)*V(k);elseif(nr(L)== n & mline(L)==1), k=nl(L);YV = YV + Ybus(n,k)*V(k);end end
Sc = conj(V(n))*(Ybus(n,n)*V(n)+ YV);Sc = conj(Sc);DP(n)= P(n)imag(Sc);if kb(n)== 1 S(n)=Sc;P(n)= real(Sc);Q(n)= imag(Sc);DP(n)=0;DQ(n)=0;Vc(n)= V(n);elseif kb(n)== 2 Q(n)= imag(Sc);S(n)= P(n)+ j*Q(n);
if Qmax(n)~= 0 Qgc = Q(n)*basemva + Qd(n)0.005;DV(n)=DV(n)+.005;end else, end else,end else,end end
if kb(n)~= 1 Vc(n)=(conj(S(n))/conj(V(n))VcI^2);Vc(n)= VcR + j*VcI;V(n)= V(n)+ accel*(Vc(n)-V(n));end end
maxerror=max(max(abs(real(DP))), max(abs(imag(DQ))));if iter == maxiter & maxerror > accuracy fprintf('nWARNING: Iterative solution did not converged after ')fprintf('%g', iter), fprintf(' iterations.nn')fprintf('Press Enter to terminate the iterations and print the results n')converge = 0;pause, else, end
end
if converge ~= 1 tech=(' ITERATIVE SOLUTION DID NOT CONVERGE');else, tech=(' Power Flow Solution by Gauss-Seidel Method');end k=0;for n = 1:nbus Vm(n)= abs(V(n));deltad(n)= angle(V(n))*180/pi;if kb(n)== 1 S(n)=P(n)+j*Q(n);Pg(n)= P(n)*basemva + Pd(n);Qg(n)= Q(n)*basemva + Qd(n)Qsh(n);end
yload(n)=(Pd(n)-j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);end
Pgt = sum(Pg);Qgt = sum(Qg);Pdt = sum(Pd);Qdt = sum(Qd);Qsht = sum(Qsh);busdata(:,3)=Vm';busdata(:,4)=deltad';clear AcurBusDPDQDVLScVcVcIVcRYVconvergedelta
2.lfybus.m程序清單
% This program obtains the Bus Admittance Matrix for power flow solution j=sqrt(-1);i = sqrt(-1);nl = linedata(:,1);nr = linedata(:,2);R = linedata(:,3);X = linedata(:,4);Bc = j*linedata(:,5);a = linedata(:, 6);nbr=length(linedata(:,1));nbus = max(max(nl), max(nr));Z = R + j*X;y= ones(nbr,1)./Z;%支路導(dǎo)納 for n = 1:nbr if a(n)<= 0 a(n)= 1;elseend
Ybus=zeros(nbus,nbus);% 將Ybus初始化為0 %非對角元素的數(shù)值
Ybus(nl(k),nr(k))=Ybus(nl(k),nr(k))-y(k)/a(k);Ybus(nr(k),nl(k))=Ybus(nl(k),nr(k));end end
% 對角元素的數(shù)值 for n=1:nbus for k=1:nbr if nl(k)==n Ybus(n,n)= Ybus(n,n)+y(k)/(a(k)^2)+ Bc(k);elseif nr(k)==n Ybus(n,n)= Ybus(n,n)+y(k)+Bc(k);else, end end end
clear Pgg
3.busout.m程序清單
% This program prints the power flow solution in a tabulated form % on the screen.disp(tech)fprintf(' Maximum Power Mismatch = %g n', maxerror)fprintf(' No.of Iterations = %g nn', iter)head =[' Bus Voltage Angle------Load---------Generation---Injected' ' No.Mag.Degree MW Mvar MW Mvar Mvar ' ' '];disp(head)for n=1:nbus fprintf(' %5g', n), fprintf(' %7.3f', Vm(n)), fprintf(' %8.3f', deltad(n)), fprintf(' %9.3f', Pd(n)), fprintf(' %9.3f', Qd(n)), fprintf(' %9.3f', Pg(n)), fprintf(' %9.3f ', Qg(n)), fprintf(' %8.3fn', Qsh(n))end
fprintf(' n'), fprintf(' Total ')fprintf(' %9.3f', Pdt), fprintf(' %9.3f', Qdt), fprintf(' %9.3f', Pgt), fprintf(' %9.3f', Qgt), fprintf(' %9.3fnn', Qsht)
4.lineflow.m程序清單
% This program is used in conjunction with lfgauss or lfNewton % for the computation of line flow and line losses.SLT = 0;fprintf('n')fprintf(' Line Flow and Losses nn')fprintf('--Line--Power at bus & line flow--Line loss--Transformern')fprintf('from to MW Mvar MVA MW Mvar tapn')for n = 1:nbus busprt = 0;for L = 1:nbr;if busprt == 0 fprintf(' n'), fprintf('%6g', n), fprintf(' %9.3f', P(n)*basemva)fprintf('%9.3f', Q(n)*basemva), fprintf('%9.3fn', abs(S(n)*basemva))
busprt = 1;else, end
if nl(L)==n k = nr(L);In =(V(n)V(n)/a(L))*y(L)+ Bc(L)*V(k);Snk = V(n)*conj(In)*basemva;Skn = V(k)*conj(Ik)*basemva;SL = Snk + Skn;SLT = SLT + SL;elseif nr(L)==n k = nl(L);In =(V(n)a(L)*V(n))*y(L)/a(L)^2 + Bc(L)/a(L)^2*V(k);Snk = V(n)*conj(In)*basemva;Skn = V(k)*conj(Ik)*basemva;SL = Snk + Skn;SLT = SLT + SL;else, end
if nl(L)==n | nr(L)==n fprintf('%12g', k), fprintf('%9.3f', real(Snk)), fprintf('%9.3f', imag(Snk))fprintf('%9.3f', abs(Snk)), fprintf('%9.3f', real(SL)), if nl(L)==n & a(L)~= 1 fprintf('%9.3f', imag(SL)), fprintf('%9.3fn', a(L))else, fprintf('%9.3fn', imag(SL))end
else, end end end
SLT = SLT/2;fprintf(' n'), fprintf(' Total loss ')fprintf('%9.3f', real(SLT)), fprintf('%9.3fn', imag(SLT))clear IkInSLSLTSknSnk
5.lfnewton.m程序清單
%Power flow solution by Newton-Raphson method ns=0;ng=0;Vm=0;delta=0;yload=0;deltad=0;nbus = length(busdata(:,1));kb=[];Vm=[];delta=[];Pd=[];Qd=[];Pg=[];Qg=[];Qmin=[];Qmax=[];Pk=[];P=[];Qk=[];Q=[];S=[];V=[];for k=1:nbus n=busdata(k,1);kb(n)=busdata(k,2);Vm(n)=busdata(k,3);delta(n)=busdata(k, 4);Pd(n)=busdata(k,5);Qd(n)=busdata(k,6);Pg(n)=busdata(k,7);Qg(n)= busdata(k,8);Qmin(n)=busdata(k, 9);Qmax(n)=busdata(k, 10);Qsh(n)=busdata(k, 11);if Vm(n)<= 0 Vm(n)= 1.0;V(n)= 1 + j*0;else delta(n)= pi/180*delta(n);V(n)= Vm(n)*(cos(delta(n))+ j*sin(delta(n)));P(n)=(Pg(n)-Pd(n))/basemva;Q(n)=(Qg(n)-Qd(n)+ Qsh(n))/basemva;S(n)= P(n)+ j*Q(n);end end
for k=1:nbus if kb(k)== 1, ns = ns+1;else, end if kb(k)== 2 ng = ng+1;else, end ngs(k)= ng;nss(k)= ns;end
Ym=abs(Ybus);t = angle(Ybus);m=2*nbus-ng-2*ns;maxerror = 1;converge=1;iter = 0;
mline=ones(nbr,1);for k=1:nbr for m=k+1:nbr if((nl(k)==nl(m))&(nr(k)==nr(m)));mline(m)=2;elseif((nl(k)==nr(m))&(nr(k)==nl(m)));mline(m)=2;else, end end end
%雅可比矩陣 clear ADCJDX
while maxerror >= accuracy & iter <= maxiter for ii=1:m for k=1:m A(ii,k)=0;%初始化雅可比矩陣 end, end
iter = iter+1;for n=1:nbus nn=n-nss(n);lm=nbus+n-ngs(n)-nss(n)-ns;J11=0;J22=0;J33=0;J44=0;for ii=1:nbr if mline(ii)==1 if nl(ii)== n | nr(ii)== n if nl(ii)== n , l = nr(ii);end if nr(ii)== n , l = nl(ii);end
J11=J11+ Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)-delta(n)+ delta(l));J33=J33+ Vm(n)*Vm(l)*Ym(n,l)*cos(t(n,l)-delta(n)+ delta(l));if kb(n)~=1 J22=J22+ Vm(l)*Ym(n,l)*cos(t(n,l)-delta(n)+ delta(l));J44=J44+ Vm(l)*Ym(n,l)*sin(t(n,l)-delta(n)+ delta(l));else, end
if kb(n)~= 1 & kb(l)~=1 lk = nbus+l-ngs(l)-nss(l)-ns;ll = l-nss(l);% J1的非對角元素
A(nn, ll)=-Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)-delta(n)+ delta(l));if kb(l)== 0 % J2的非對角元素 A(nn, lk)=Vm(n)*Ym(n,l)*cos(t(n,l)-delta(n)+ delta(l));end if kb(n)== 0 % J3的非對角元素
A(lm, ll)=-Vm(n)*Vm(l)*Ym(n,l)*cos(t(n,l)-delta(n)+delta(l));end
if kb(n)== 0 & kb(l)== 0 % J4的非對角元素
A(lm, lk)=-Vm(n)*Ym(n,l)*sin(t(n,l)-delta(n)+ delta(l));end elseend else , end else, end end
Pk = Vm(n)^2*Ym(n,n)*cos(t(n,n))+J33;Qk =-Vm(n)^2*Ym(n,n)*sin(t(n,n))-J11;if kb(n)== 1 P(n)=Pk;Q(n)= Qk;end% Swing bus P if kb(n)== 2 Q(n)=Qk;if Qmax(n)~= 0 Qgc = Q(n)*basemva + Qd(n)0.01;end else, end else,end else,end end
if kb(n)~= 1 A(nn,nn)= J11;% J1對角元素 DC(nn)= P(n)-Pk;end
if kb(n)== 0 A(nn,lm)= 2*Vm(n)*Ym(n,n)*cos(t(n,n))+J22;% J2對角元素 A(lm,nn)= J33;% J3對角元素
A(lm,lm)=-2*Vm(n)*Ym(n,n)*sin(t(n,n))-J44;% J4對角元素 DC(lm)= Q(n)-Qk;end end
DX=ADC';for n=1:nbus nn=n-nss(n);lm=nbus+n-ngs(n)-nss(n)-ns;if kb(n)~= 1 delta(n)= delta(n)+DX(nn);end if kb(n)== 0 Vm(n)=Vm(n)+DX(lm);end end
maxerror=max(abs(DC));if iter == maxiter & maxerror > accuracy fprintf('nWARNING: Iterative solution did not converged after ')fprintf('%g', iter), fprintf(' iterations.nn')fprintf('Press Enter to terminate the iterations and print the results n')converge = 0;pause, else, end
end
if converge ~= 1 tech=(' ITERATIVE SOLUTION DID NOT CONVERGE');else, tech=(' Power Flow Solution by Newton-Raphson Method');end
V = Vm.*cos(delta)+j*Vm.*sin(delta);deltad=180/pi*delta;i=sqrt(-1);k=0;for n = 1:nbus if kb(n)== 1 k=k+1;S(n)= P(n)+j*Q(n);Pg(n)= P(n)*basemva + Pd(n);Qg(n)= Q(n)*basemva + Qd(n)Qsh(n);Pgg(k)=Pg(n);Qgg(k)=Qg(n);end
yload(n)=(Pd(n)-j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);end
busdata(:,3)=Vm';busdata(:,4)=deltad';Pgt = sum(Pg);Qgt = sum(Qg);Pdt = sum(Pd);Qdt = sum(Qd);Qsht = sum(Qsh);
6.decouple.m程序清單
% Fast Decoupled Power Flow Solution ns=0;Vm=0;delta=0;yload=0;deltad=0;nbus = length(busdata(:,1));kb=[];Vm=[];delta=[];Pd=[];Qd=[];Pg=[];Qg=[];Qmin=[];Qmax=[];Pk=[];P=[];Qk=[];Q=[];S=[];V=[];for k=1:nbus n=busdata(k,1);kb(n)=busdata(k,2);Vm(n)=busdata(k,3);delta(n)=busdata(k, 4);Pd(n)=busdata(k,5);Qd(n)=busdata(k,6);Pg(n)=busdata(k,7);Qg(n)= busdata(k,8);Qmin(n)=busdata(k, 9);Qmax(n)=busdata(k, 10);Qsh(n)=busdata(k, 11);if Vm(n)<= 0 Vm(n)= 1.0;V(n)= 1 + j*0;else delta(n)= pi/180*delta(n);V(n)= Vm(n)*(cos(delta(n))+ j*sin(delta(n)));P(n)=(Pg(n)-Pd(n))/basemva;Q(n)=(Qg(n)-Qd(n)+ Qsh(n))/basemva;S(n)= P(n)+ j*Q(n);end
if kb(n)== 1, ns = ns+1;else, end nss(n)= ns;end
Ym = abs(Ybus);t = angle(Ybus);ii=0;for ib=1:nbus if kb(ib)== 0 | kb(ib)== 2 ii = ii+1;jj=0;for jb=1:nbus if kb(jb)== 0 | kb(jb)== 2 jj = jj+1;B1(ii,jj)=imag(Ybus(ib,jb));else,end end
else, end end
ii=0;for ib=1:nbus if kb(ib)== 0 ii = ii+1;jj=0;for jb=1:nbus if kb(jb)== 0 jj = jj+1;B2(ii,jj)=imag(Ybus(ib,jb));else,end end
else, end end
B1inv=inv(B1);B2inv = inv(B2);
maxerror = 1;converge = 1;iter = 0;
mline=ones(nbr,1);for k=1:nbr for m=k+1:nbr if((nl(k)==nl(m))&(nr(k)==nr(m)));mline(m)=2;elseif((nl(k)==nr(m))&(nr(k)==nl(m)));mline(m)=2;else, end end end
% 開始迭代
while maxerror >= accuracy & iter <= maxiter % 檢驗不平衡功率 iter = iter+1;id=0;iv=0;for n=1:nbus nn=n-nss(n);J11=0;J33=0;for ii=1:nbr if mline(ii)==1 if nl(ii)== n | nr(ii)== n if nl(ii)== n, l = nr(ii);end if nr(ii)== n, l = nl(ii);end
J11=J11+ Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)-delta(n)+ delta(l));J33=J33+ Vm(n)*Vm(l)*Ym(n,l)*cos(t(n,l)-delta(n)+ delta(l));else , end else, end end
Pk = Vm(n)^2*Ym(n,n)*cos(t(n,n))+J33;Qk =-Vm(n)^2*Ym(n,n)*sin(t(n,n))-J11;if kb(n)== 1 P(n)=Pk;Q(n)= Qk;end% Swing bus P if kb(n)== 2 Q(n)=Qk;Qgc = Q(n)*basemva + Qd(n)0.005;end% the specified limits.else, end else,end else,end end
if kb(n)~= 1 id = id+1;DP(id)= P(n)-Pk;DPV(id)=(P(n)-Pk)/Vm(n);end
if kb(n)== 0 iv=iv+1;DQ(iv)= Q(n)-Qk;DQV(iv)=(Q(n)-Qk)/Vm(n);end end
Dd=-B1DPV';DV=-B2DQV';id=0;iv=0;for n=1:nbus if kb(n)~= 1 id = id+1;delta(n)= delta(n)+Dd(id);end if kb(n)== 0 iv = iv+1;Vm(n)=Vm(n)+DV(iv);end end
maxerror=max(max(abs(DP)),max(abs(DQ)));if iter == maxiter & maxerror > accuracy fprintf('nWARNING: Iterative solution did not converged after ')fprintf('%g', iter), fprintf(' iterations.nn')fprintf('Press Enter to terminate the iterations and print the results n')converge = 0;pause, else, end
end
if converge ~= 1 tech=(' ITERATIVE SOLUTION DID NOT CONVERGE');else, tech=(' Power Flow Solution by Fast Decoupled Method');end k=0;V = Vm.*cos(delta)+j*Vm.*sin(delta);deltad=180/pi*delta;clear A;clear DC;clear DX i=sqrt(-1);for n = 1:nbus if kb(n)== 1 S(n)=P(n)+j*Q(n);Pg(n)= P(n)*basemva + Pd(n);Qg(n)= Q(n)*basemva + Qd(n)Qsh(n);k=k+1;Pgg(k)=Pg(n);end
yload(n)=(Pd(n)-j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);end
busdata(:,3)=Vm';busdata(:,4)=deltad';Pgt = sum(Pg);Qgt = sum(Qg);Pdt = sum(Pd);Qdt = sum(Qd);Qsht = sum(Qsh);clear PkQkDPDQJ11J33B1B1invB2B2invDPVDQVDddeltaibidiiivjbjj
四、30節(jié)點電力系統(tǒng)計算實例
潮流計算時，必須將前面的六個程序保存在MATLAB目錄下格式為.m的文件，然后在MATLAB的命令窗口輸入如下命令： clear basemva = 100;accuracy = 0.001;accel = 1.8;maxiter = 100;% 30節(jié)點電力系統(tǒng)
% 母線--母線 電壓 相角 負(fù)載 發(fā)電機 注入功率
% 編號節(jié)點 幅值 角度有功 無功有功 無功 無功最小值 無功最大值 無功 busdata=[1 1 1.06 0.0 0.0 0.0 0.0 0.0 0 0 0 2 2 1.043 0.0 21.70 12.7 40.0 0.0-40 50 0 3 0 1.0 0.0 2.4 1.2 0.0 0.0 0 0 0 4 0 1.06 0.0 7.6 1.6 0.0 0.0 0 0 0 5 2 1.01 0.0 94.2 19.0 0.0 0.0-40 40 0 6 0 1.0 0.0 0.0 0.0 0.0 0.0 0 0 0 7 0 1.0 0.0 22.8 10.9 0.0 0.0 0 0 0 8 2 1.01 0.0 30.0 30.0 0.0 0.0-30 40 0 9 0 1.0 0.0 0.0 0.0 0.0 0.0 0 0 0 10 0 1.0 0.0 5.8 2.0 0.0 0.0-6 24 19 11 2 1.082 0.0 0.0 0.0 0.0 0.0 0 0 0 12 0 1.0 0 11.2 7.5 0 0 0 0 0 13 2 1.071 0 0 0.0 0 0-6 24 0 14 0 1 0 6.2 1.6 0 0 0 0 0 15 0 1 0 8.2 2.5 0 0 0 0 0 16 0 1 0 3.5 1.8 0 0 0 0 0 17 0 1 0 9.0 5.8 0 0 0 0 0 18 0 1 0 3.2 0.9 0 0 0 0 0 19 0 1 0 9.5 3.4 0 0 0 0 0 20 0 1 0 2.2 0.7 0 0 0 0 0 21 0 1 0 17.5 11.2 0 0 0 0 0 22 0 1 0 0 0.0 0 0 0 0 0 23 0 1 0 3.2 1.6 0 0 0 0 0 24 0 1 0 8.7 6.7 0 0 0 0 4.3 25 0 1 0 0 0.0 0 0 0 0 0 26 0 1 0 3.5 2.3 0 0 0 0 0 27 0 1 0 0 0.0 0 0 0 0 0 28 0 1 0 0 0.0 0 0 0 0 0 29 0 1 0 2.4 0.9 0 0 0 0 0 30 0 1 0 10.6 1.9 0 0 0 0 0];% 線路數(shù)據(jù)
% bus bus R X 1/2 B 1 for lines linedata=[1 2 0.0192 0.0575 0.02640 1 1 3 0.0452 0.1852 0.02040 1 2 4 0.0570 0.1737 0.01840 1 3 4 0.0132 0.0379 0.00420 1 2 5 0.0472 0.1983 0.02090 1 2 6 0.0581 0.1763 0.01870 1 4 6 0.0119 0.0414 0.00450 1 5 7 0.0460 0.1160 0.01020 1 6 7 0.0267 0.0820 0.00850 1 6 8 0.0120 0.0420 0.00450 1 6 9 0.0 0.2080 0.0 0.978 6 10 0.5560 0 0.969 9 11 0.2080 0 1 9 10 0.1100 0 1 4 12 0.2560 0 0.932 12 13 0.1400 0 1 12 14.1231.2559 0 1 12 15.0662.1304 0 1 12 16.0945.1987 0 1 14 15.2210.1997 0 1 16 17.0824.1923 0 1 15 18.1073.2185 0 1 18 19.0639.1292 0 1 19 20.0340.0680 0 1 10 20.0936.2090 0 1 10 17.0324.0845 0 1 10 21.0348.0749 0 1 10 22.0727.1499 0 1 21 22.0116.0236 0 1 15 23.1000.2020 0 1 22 24.1150.1790 0 1 23 24.1320.2700 0 1 24 25.1885.3292 0 1 25 26.2544.3800 0 1 25 27.1093.2087 0 1 28 27 0.3960 0 0.968 27 29.2198.4153 0 1 27 30.3202.6027 0 1 29 30.2399.4533 0 1 8 28.0636.2000 0.0214 1 6 28.0169.0599 0.065 1];最后運行程序輸入以下命令： lfybus % 形成節(jié)點導(dǎo)納矩陣 lfgauss % 高斯-賽德爾法潮流計算 busout % 屏幕顯示潮流計算結(jié)果 lineflow % 計算并顯示線路潮流和損耗
將lfgauss變?yōu)閘fnewton/decouple,即可使用牛頓-拉夫遜法/快速解耦法進行潮流計算，輸入以上4個命令行后，即可得到潮流計算結(jié)果：
第二篇：電力系統(tǒng)潮流計算程序設(shè)計
電力系統(tǒng)潮流計算程序設(shè)計
姓名：韋應(yīng)順
學(xué)號：2011021052 電力工程學(xué)院
牛頓—拉夫遜潮流計算方法具有能夠?qū)⒎蔷€性方程線性化的特點，而使用MATLAB語言是由于MATLAB語言的數(shù)學(xué)邏輯強，易編譯。
【】【】1.MATLAB程序12
Function tisco %這是一個電力系統(tǒng)潮流計算的程序 n=input（‘n請輸入節(jié)點數(shù)：n=’）； m=input(‘請輸入支路數(shù)：m=’);ph=input(‘n請輸入平衡母線的節(jié)點號：ph=’)； B1=input（‘n請輸入支路信號：B1=’）;%它以矩陣形式存貯支路的情況，每行存貯一條支路 %第一列存貯支路的一個端點 %第二列存貯支路的另一個端點 %第三列存貯支路阻抗
%第四列存貯支路的對地導(dǎo)納
%第五列存貯變壓器的變比，注意支路為1 %第六列存貯支路的序號
B2=input(‘n請輸入節(jié)點信息：B2=’)； %第一列為電源側(cè)的功率 %第二列為負(fù)荷側(cè)的功率 %第三列為該點的電壓值
%第四列為該點的類型：1為PQ，2為PV節(jié)點，3為平衡節(jié)點 A=input（‘n請輸入節(jié)點號及對地阻抗：A=’）； ip=input（‘n請輸入修正值：ip=’）； %ip為修正值);Y=zeros（n）；
Y(p,q)=Y(p,q)-1./(B1(i3)*B1(i5);e=zeros（1，n）；
Y(p,q)=Y(p,q);f=zeros（1，n）；
no=2*ph=1； Y(q,q)=Y(q,q)+1./B1(i3)+B1(i4)/2；
End for i=1：n
G=real(Y);if A（i2）=0
B=imag(Y);p=A(i1)；
Y(p p)=1./A(i2);for i=1：n End e(i)=real(B2(i3));End f(i)=imag(B2(i3));For i=1:m S(i)=B2(i1)-B2(i2);p=B1(i1);V(i)=B2(i3);p=B1(i2);end Y(p,p)=Y(p,p)+1./(B1(i3)*B1(i5)^2+B1(i4)./2P=real(S);Q=imag(S);[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no);J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no);[De,Di]=hxf(J,D,F,ph,n,no);t=0;while
max(abs(De))>ip&max(abs(Dfi)>ip
t=t+1;
e=e+De；
f=f+Df;
[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)；
J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no)；
[De,Df]=hxf(J,Df,ph,n,no);end v=e+f*j;for i=1:n hh(i)=conj(Y(ph,i)*v(i));end S(ph)=sum(hh)*v(ph);B2(ph,1)=S(ph);V=abs(v);
jd=angle(v)*180/p;resulte1=[A(:,1),real(v),imag(v),V,jd,real(S’)，imag(S’),real(B2(:1)),imag(B2(:1)),real(B2(:2)),imag(B2(:,2))];for i=1:m
a(i)=conj((v(B1(i1))/B1(i5)-v(B1(i2))/B1(i3));
b(i)=v(B1(i1))*a(i)-j*B1(i4)*v(B1(i))^2/2;
c(i)=-v(B1(i2))*a(i)-j*B1(i4)*v(B1(i2))^2/2;end result2=[B1(:,6),B1(:,1),B1(:,2),real(b’),imag(b’),real(c’),imag(c’), real(b’+c’),imag(b’+c’)];printcut(result1,S,b,c,result2);type resultm function [C,D,Df]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)%該子程序是用來求取Df for i=1:n
If
i=ph
C(i)=0;
D(i)=0;
For j=i:n
C(i)=C(i)+G(i,j)*e(j)-B(i,j)*f(j);D(i)=D(i)+G(i,j)*f(j)+B(i,j)*e(j);end
P1=C(i)*e(i)+D(i)*f(i);Q1=C(i)*f(i)-D(i)*e(i);V1=e(i)^2+f(i)^2;If
B2(i4)=2 p=2*i-1;
Df(p)=P(i)-P1;p=p+1;else p=2*i-1;
Df(p)=P(i)-P1;p=p+1;
Df(p)=Q(i)-Q1;end end end Df=Df’;If ph=n Df(no?=[];end
function [De,Df]=hxf(J,Df,ph,n,no)%該子函數(shù)是為求取De Df DX=JDf;DX1=DX;
x1=length(DX1);if ph=n DX(no)=0;DX(no+1)=0;
For i=(no+2):(x1+2)DX(i)=DX1(i-2);End Else
DX=[DX1,0,0];End k=0;
[x,y]=size(DX);For i=1:2:x K=k+1;
Df(k)=DX(i);De(k)=DX(i+1);End End case 2 Function for j=1:n J=jacci(Y,G,B,PQ,e,f,V,C,D,B2,n,ph,no)X1=G(i,j)*f(i)-B(i,j)*e(i);
X2=G(i,j)*e(i)+B(i,j)*f(i);%該子程序是用來求取jacci矩陣
for i=1:n X3=0;switch B2(i4)X4=0;case 3 P=2*i-1;continue q=2*j-1;case 1 J(p,q)=X1;for j=1:n m=p+1;if
J=&J=ph J(m,q)=X3;X1=G(i)*f(i)-B(i,j)*e(i);q=q+1;X2=G(i,j)*e(i)+B(i,j)*f(i);J(p,q)=X2;X3=-X2;J(m,q)=X4;X4=X1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);p=2*i-1;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=2*j-1;X3=0;J(p,q)=X1;X4=0;m=p+1;P=2*i-1;J(p,q)=X2;q=2*j-1;J(m,q)=X4;J(p,q)=X1;Else if j=&j=jph m=p+1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X3;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=q+1;X3= C(i)+G(i,j)*e(i)-B(i,j)*f(i);J(p,q)=X2;X4= C(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X4;P=2*i-1;end q=2*j-1;end J(p,q)=X1;end m=p+1;end J(m,q)=X3;if ph=n q=q+1;J(no:)=[];J(p,q)=X2;J(no:)=[];J(m,q)=X4;J(:,no)=[];End J(:,no)=[];End
2實例驗證 【例題】設(shè)有一系統(tǒng)網(wǎng)絡(luò)結(jié)線見圖1，各支路阻抗和各節(jié)點功率均已以標(biāo)幺值標(biāo)示于圖1中，其中節(jié)點2連接的是發(fā)電廠，設(shè)節(jié)點1電壓保持U1＝1.06定值，試計算其中的潮流分布，請輸入節(jié)點數(shù)：n=5 請輸入支路數(shù)：m=7 請輸入平衡母線的節(jié)點號：ph=l 請輸入支路信息：
BI=[ l 2 0．02+0．06i O l 1；1 3 0．08+0．24i 0 1 2；2 3 0．06+0．18i 0 l 3： 2 4 0．06+0．18i O l 4： 2 5 0．04+0．12i 0 l 5： 3 4 0．01+0．03i 0 l 6： 4 5 0．08+0．24i O 1 7] 請輸入節(jié)點信息：
B2=[ 0 0 1．06 3；0．2+0．20i 0 1 1；一O．45一O．15i 0 l l；一0．4-0．05i 0 l 1；一0．6—0．1i 0 1 l] 請輸入節(jié)點號及對地阻抗： A=[l 0；2 0；3 0；4 0；5 O ] 請輸入修正值：ip=0．000 0l
參考文獻
[1]陳珩．電力系統(tǒng)穩(wěn)定分析[M]．北京：中國電力出版社，2002：139—187．
[2]鄭阿奇．MATLAB實用教程[M]．北京：電子工業(yè)出版社，2005：1-243．
[3] 束洪春，孫士云，等．云電送粵交商流混聯(lián)系統(tǒng)全過 程動態(tài)電壓研究[J】．中國電力，2008，4l(10)：l-4． SHU Hong—ch吼，SUN Shi-yun，et a1．Research on fun prc'cess dyn鋤ic Voltage stabil時of hybrid AC／DC poWer tmnsmission System舶m Yu衄an proVince to G啪gdong province【J】．Electric Power，2008，4l(10)： l-4．
[4] 朱新立，湯涌，等．大電網(wǎng)安全分析的全過程動態(tài)仿 真技術(shù)[J】．電網(wǎng)技術(shù)，2008，32(22)：23—28． SONG Xin—Ii，TANG Yof唱，et a1． Full dyn鋤ic simulation for the stabilhy a眥lysis of large power system【J】．Power System融IlrIolo影，2008，32(22)： 23．28．
[5]Roytelm鋤I，Shallidehpour S M．A comprehcnsivc long teml dynaIIlic simulation for powcr system recoVery【J】． IEEE Transactions 0n Power Systems，1994，9(3)． [6] 石雩梅，汪志宏，等．發(fā)電機勵磁系統(tǒng)數(shù)學(xué)模型及參 數(shù)對電網(wǎng)動態(tài)穩(wěn)定性分析結(jié)果影響的研究[J】．繼電 器，2007，35(21)：22-27．
SHI Xue．mei，WANG Zlli-hon舀et a1．Iksearch on the innuence of g鋤e翰to璐baScd ∞de詛iled excitation system models柚d parameterS t0 power鏟id dyn鋤ic stabil時【J】．Relay，2007，35(2 1)：22-27．
[7] 方思立，朱方．快速勵磁系統(tǒng)對系統(tǒng)穩(wěn)定的影響[J】．中 國電機工程學(xué)報，1986，6(1)：20．28．
FANG Si．1i，ZHU Fang．The effbct of f弧t．respon∞
excitation system on the stability of power netwofk【J】． Proceedings ofthe CSEE，1986，6(1)：20-28．
[8] 劉取．電力系統(tǒng)穩(wěn)定性及發(fā)電機勵磁控制[M】．北京： 中國電力出版社，2007．
LIU Qu．Power system S詛bility鋤d generator excitation control【M】．BeUing：ChiIla Electric Powef Press，2007． [9] Dallachy J L，Anderson T．EXperience with rcplacing ro詛ting exciters wim static exciters【J】．1k InStitution of Electrical Engineers，1 996．
[10] 陳利芳，陳天祿．淺談自并勵勵磁系統(tǒng)在大容量機組 中的應(yīng)用【J】．繼電器，2007，35(1)：8l培4． CHEN Li-f抽島CHEN Tian—lIL Application of 辯l仁exci組tion mode in large capacity髫memtor unit【J】． ReIay'2007，35(1)：81-84．
[11] 方思立，劉增煌，孟慶和．大型汽輪發(fā)電機自并勵勵 磁系統(tǒng)的應(yīng)用條件【J]．中國電力，1994，27(12)：61．63． FANG Si．Ii，LIU Zeng-hu鋤g，MENG Qin爭hc．m application conditions of large turbine generator self-excitation system【J】．Electric Powef，1994，27(12)： 61．63．
[12]梁小冰，黃方能．利用EMTDC進行長持續(xù)時間過程 的仿真研究【J】．電網(wǎng)技術(shù)，2002，26(9)：55．57． LIANG Xiao-bing，HUANG Fan爭眥ng．How to cany out simulalion of long dul‘a(chǎn)tion processes by use of EMTDC【J】．Power System 11echnology，2002，26(9)： 55-57．
[13]王卉，陳楷，彭哲，等．?dāng)?shù)字仿真技術(shù)在電力系統(tǒng)中 的應(yīng)用及常用的幾種數(shù)字仿真工具【J】．繼電器，2004，32(21)：7l一75．
wANG Hui，CHEN Kai，PENG zhe，et a1．Application of digital simulation眥hniques棚d severaJ simulation tools in power system[J】．Relay，2004，32(21)：71·75．
[14]IEEE Power Engmeering Socie哆．IEEE std 421．5．2005 IEEE玎ccOmmended practice for excitation system models for power system stabiI時studies【s】．
第三篇：用matlab電力系統(tǒng)潮流計算
題目：潮流計算與matlab
教學(xué)單位 電氣信息學(xué)院
姓 名 學(xué) 號
年 級
專 業(yè) 電氣工程及其自動化
指導(dǎo)教師
職 稱 副教授
摘 要
電力系統(tǒng)穩(wěn)態(tài)分析包括潮流計算和靜態(tài)安全分析。本文主要運用的事潮流計算，潮流計算是電力網(wǎng)絡(luò)設(shè)計與運行中最基本的運算，對電力網(wǎng)絡(luò)的各種設(shè)計方案及各種運行方式進行潮流計算，可以得到各種電網(wǎng)各節(jié)點的電壓，并求得網(wǎng)絡(luò)的潮流及網(wǎng)絡(luò)中的各元件的電力損耗，進而求得電能損耗。本位就是運用潮流計算具體分析，并有MATLAB仿真。
關(guān)鍵詞： 電力系統(tǒng) 潮流計算 MATLAB
Abstract Electric power system steady flow calculation and analysis of the static safety analysis.This paper, by means of the calculation, flow calculation is the trend of the power network design and operation of the most basic operations of electric power network, various design scheme and the operation ways to tide computation, can get all kinds of each node of the power grid voltage and seek the trend of the network and the network of the components of the power loss, and getting electric power.The standard is to use the power flow calculation and analysis, the specific have MATLAB simulation.Key words: Power system;Flow calculation;MATLAB simulation
目 錄 任務(wù)提出與方案論證....................................................................................................................................2 2 總體設(shè)計........................................................................................................................................................3 2.1潮流計算等值電路.............................................................................................................................3 2.2建立電力系統(tǒng)模型.............................................................................................................................3 2.3模型的調(diào)試與運行.............................................................................................................................3 3 詳細(xì)設(shè)計........................................................................................................................................................4 3.1 計算前提............................................................................................................................................4 3.2手工計算.............................................................................................................................................7 4設(shè)計圖及源程序...........................................................................................................................................11 4.1MATLAB仿真.......................................................................................................................................11 4.2潮流計算源程序...............................................................................................................................11 5 總結(jié).............................................................................................................................................................31 參考文獻..........................................................................................................................................................32 任務(wù)提出與方案論證
潮流計算是在給定電力系統(tǒng)網(wǎng)絡(luò)結(jié)構(gòu)、參數(shù)和決定系統(tǒng)運行狀態(tài)的邊界條件的情況下確定系統(tǒng)穩(wěn)態(tài)運行狀態(tài)的一種基本方法,是電力系統(tǒng)規(guī)劃和運營中不可缺少的一個重要組成部分。可以說,它是電力系統(tǒng)分析中最基本、最重要的計算,是系統(tǒng)安全、經(jīng)濟分析和實時控制與調(diào)度的基礎(chǔ)。常規(guī)潮流計算的任務(wù)是根據(jù)給定的運行條件和網(wǎng)路結(jié)構(gòu)確定整個系統(tǒng)的運行狀態(tài)，如各母線上的電壓（幅值及相角）、網(wǎng)絡(luò)中的功率分布以及功率損耗等。潮流計算的結(jié)果是電力系統(tǒng)穩(wěn)定計算和故障分析的基礎(chǔ)。在電力系統(tǒng)運行方式和規(guī)劃方案的研究中，都需要進行潮流計算以比較運行方式或規(guī)劃供電方案的可行性、可靠性和經(jīng)濟性。同時，為了實時監(jiān)控電力系統(tǒng)的運行狀態(tài)，也需要進行大量而快速的潮流計算。因此，潮流計算是電力系統(tǒng)中應(yīng)用最廣泛、最基本和最重要的一種電氣運算。在系統(tǒng)規(guī)劃設(shè)計和安排系統(tǒng)的運行方式時，采用離線潮流計算；在電力系統(tǒng)運行狀態(tài)的實時監(jiān)控中，則采用在線潮流計算。是電力系統(tǒng)研究人員長期研究的一個課題。它既是對電力系統(tǒng)規(guī)劃設(shè)計和運行方式的合理性、可靠性及經(jīng)濟性進行定量分析的依據(jù)，又是電力系統(tǒng)靜態(tài)和暫態(tài)穩(wěn)定計算的基礎(chǔ)。
潮流計算經(jīng)歷了一個由手工到應(yīng)用數(shù)字電子計算機的發(fā)展過程，現(xiàn)在的潮流算法都以計算機的應(yīng)用為前提用計算機進行潮流計算主要步驟在于編制計算機程序，這是一項非常復(fù)雜的工作。對系統(tǒng)進行潮流分析，本文利用 MATLAB中的SimpowerSystems工具箱設(shè)計電力系統(tǒng)，在simulink 環(huán)境下，不僅可以仿真系統(tǒng)的動態(tài)過程，還可以對系統(tǒng)進行穩(wěn)態(tài)潮流分析。
總體設(shè)計
SimpowerSystems使用Simulink環(huán)境，可以將該系統(tǒng)中的發(fā)電機、變壓器，線路等模型聯(lián)結(jié)起來，形成電力系統(tǒng)仿真模擬圖。在加人測量模塊，并對各元件的參數(shù)進行設(shè)置后，用measurement和sink中的儀器可以觀察各元件的電壓、電流、功率的大小。
2.1潮流計算等值電路
10MWYN,d114?63MWVA15MWGGGG120MW10kV??p0?15.7kW??ps?73kW?I%?0.5?0?Vs%?10.5YN,d1116MWVA4?63MW“?xd?0.134?x2?0.161?x?0.06?0?cos?N?0.8510kV??p0?11kW??ps?50kW?I%?0.55?0?Vs%?10.5YN,d112?10MWVA35kV32km25MW110kV80km25MW110kV70km110kVYN,d112?20MWVA20MWGGG4?15MW”?xd?0.136?x2?0.16?x?0.073?0?cos?N?0.8??p0?18.6kW??ps?89kW?I%?0.530MW?0?Vs%?10.510kV??p0?15.7kW??ps?73kW?I%?0.5GG?0V%?10.5s?YN,d112?16MWVA63MWVA??p0?44kW??ps?121kW10kV?I%?0.35?0?Vs%?10.5GG35MWYN,Y,d112?10MVA10kVGGG3?12MW1?50MW“?xd?0.128?x2?0.154?x?0.054?0?cos?N?0.852?25MW”?xd?0.128?x2?0.157?x?0.0591?0?cos?N?0.8?x?0.136?x2?0.161?x?0.075?0?cos?N?0.8“d
2.2建立電力系統(tǒng)模型
在Simulink中按照電力系統(tǒng)原型選擇元件進行建模。所建立的模型和建立的方法在詳細(xì)設(shè)計中詳述。
在電力系統(tǒng)模型的建立工程中主要涉及到的是：元器件的選擇及其參數(shù)的設(shè)置；發(fā)電機選型；變壓器選擇；線路的選擇；負(fù)荷模型的選擇；母線選擇。
2.3模型的調(diào)試與運行
建立系統(tǒng)模型，并設(shè)置好參數(shù)以后，就可以在Simulink環(huán)境下進行仿真運行。運行的具體結(jié)果和分析也在詳細(xì)設(shè)計中詳述。
30km35kV0km31??p0?44kW??ps?121kW?I%?0.35?0?Vs%?10.5??p0?13.2kW??ps?63kW?I0%?0.55?V%?10.5s(1?2)80MW?Vs(2?3)%?6.55MW??Vs(1?3)%?17.53 詳細(xì)設(shè)計
3.1 計算前提
首先是發(fā)電機的參數(shù)計算，先對5個發(fā)電廠簡化為5臺發(fā)電機來計算。發(fā)電機G1：
P1?4?15?60MWQ1?60?tan(arccos0.8)?45MVar發(fā)電機G2：
P2?4?63?252MWQ2?252?tan(arccos0.85)?156MVar發(fā)電機G3：
P3?3?12?36MWQ3?36?tan(arccos0.8)?27MVar發(fā)電機G4：
P4?1?50?50MWQ4?50?tan(arccos0.85)?31MVar發(fā)電機G5：
P5?2?25?50MWQ5?50?tan(arccos0.8)?37.5MVar
其次是變電站的參數(shù)計算，我們還是對7個變電站簡化為7臺變壓器來計算。變壓器T1：
2?ps?VN73?11023RT1??10??103?3.450?232SN(16?10)2Vs%?VN10.5?1102XT1??10??10?79.406?SN16?103?S01??p0?j變壓器T2：（雙并聯(lián)）
I0%?SN?(0.0157?j0.0800)MVA 100RT2XT221?ps?VN189?11023???10???103?1.346?2322SN2(20?10)21Vs%?VN110.5?1102???10???10?31.7625? 32SN220?10?S02?2?(?p0?j變壓器T3：（四并聯(lián)）
I0%?SN)?(0.0372?j0.2000)MVA100 1?ps?VN21121?11023RT3???10???103?0.092?2324SN4(63?10)XT31Vs%?VN2110.5?1102???10???10?5.042? 4SN463?103I0%?SN)?(0.1760?j0.8820)MVA100?S03?4?(?p0?j變壓器T4：（雙并聯(lián)）
1RT1?1.7250?21 XT4?XT1?39.7030?2?S04?2?S01?(0.0314?j0.1600)MVART4?變壓器T5：
RT5?4RT3?0.3680?XT5?4XT3?20.168??S05?1?S03?(0.0440?j0.2205)MVA4163?3523???10?0.386? 322(10?10)
變壓器T6：（兩個三繞組變壓器并聯(lián)）
RT6?1?RT6?2?RT6?31?[Vs(1?2)%?Vs(1?3)%?Vs(2?3)%]?10.7521Vs2%??[Vs(1?2)%?Vs(2?3)%?Vs(1?3)%]??0.25
21Vs3%??[Vs(1?3)%?Vs(2?3)%?Vs(1?2)%]?6.752Vs1%?21Vs1%?VNXT6?1???10?6.584?2SNXT6?2XT6?321Vs2%?VN???10??0.153?2SN21Vs3%?VN???10?4.134? 2SN?S06?2?(?P06?j變壓器T7：（雙并聯(lián)）
I0%?10)?(0.0264?j0.1100)MVA 100 RT7XT721?ps?VN150?3523???10???103?0.306?2322SN2(10?10)21Vs%?VN110.5?352???10???10?6.431?2SN210?103
?S07?2?(?p0?jI0%?SN)?(0.0220?j0.1100)MVA100再次是傳輸線參數(shù)計算，5條傳輸線的具體計算如下。
根據(jù)教材查得r0?0.21?/km x0?0.4?/km b0?2.8?10S/km ?6線路L1：
線路L2：
線路L3：（雙回路）
線路L4：
線路L5：（雙回路）RL1?r0?l1?0.21?40?8.4?XL1?x0?l1?0.4?40?16?B?6L1?b0?l1?2.8?10?40?1.12?10?4S ?Q1L1??2BL1V2N??0.6776MVarRL2?r0?l2?0.21?130?27.3?XL2?x0?l2?0.4?130?52?B?6L2?b0?l2?2.8?10?130?3.64?10?4S ?Q1L2??2BL2V2N??2.2022MVarR?12?rl1L30?3?2?0.21?70?7.35?X11L3?2?x0?l3?2?0.4?70?14? BL3?2?b?40?l3?2?2.8?10?6?70?3.92?10S?Q1L3??2B2L3VN??2.3716MVarRL4?r0?l4?0.21?60?12.6?XL4?x0?l4?0.4?60?24?BL4?b0?l4?2.8?10?6?60?1.68?10?4S ?Q12L4??2BL4VN??1.0164MVar
11RL5??r0?l5??0.21?20?2.1?2211XL5??x0?l5??0.4?20?4? 22BL5?2?b0?l5?2?2.8?10?6?20?1.12?10?4S1?QL5??BL3VN2??0.0686MVar23.2手工計算
FLR1：
P2102?ST1?2(RT1?jXT1)?(3.450?j74.406)?(0.0285?j0.6562)MVA2VN110Sa?10MW??ST1??S01?j?QL1?(10.0442?j0.1142)MVAP2?Q210.04422?0.11422?SL1?(RL1?jXL1)?(8.4?j16)?(0.070?j0.1334)MVAVN21102?ST2P2?Q2402?452?(RT2?jXT2)?(1.346?j31.7625)?(0.4032?j9.5156)MVAVN21102FLR2Sb?SG1?20??ST2?60?j45?20?0.4032?j9.5156?(39.5968?j35.4844)MVASc?Sb?Sa?25?jQL1??SL1?(4.4826?j35.9144)MVA：
?ST3P2?Q22522?1562?(RT3?jXT3)?(0.092?j5.042)?(0.6679?j36.6024)MVA22VN110P?Q4.4931?34.1048(R?jX)?(27.3?j52)?(2.67?j5.0854)MVAL2L222VN1102222Sc'?(4.4931?j34.1048)MVA?SL2?FLRSd?SG2?Sc'?120??ST3??S03?jQL2??SL2?(132.9792?j149.229)MVA3：
?ST4P2?Q262?272?(RT4?jXT4)?(1.725?j39.703)?(0.1091?j2.5101)MVAVN21102P2?Q2133.59552?149.99562?(RL3?jXL3)?(7.35?j14)?(24.51?j46.682)MVA22VN110'Sd?(133.5955?j149.9956)MVA?SL3'Se?SG3?Sd?30?25??ST4??S04?jQL3??SL3?(89.945?j130.0151)MVAFLR4： ?ST5P2?Q2502?312?(RT5?jXT5)?(0.368?j20.168)?(0.1052?j5.7687)MVA22VN110P2?Q292.74872?133.99372?(RL4?jXL4)?(12.6?j24)?(27.654?j52.674)MVA22VN110Se'?(92.7481?j133.9937)MVA?SL4Sf?SG4?Se'?80??ST5??S05?jQL4??SL4?(34.9449?j107.3469)MVAFLR5: 152?ST7?2?(0.306?j6.431)?(0.0562?j1.1812)MVA35Sh?15??ST7??S07?j?QL5?(15.0782?j0.3422)MVA15.07822?0.34222?SL5??(2.1?j4)?(0.3899?j0.743)MVA352Si?Sh??SL5??S06?j?QL5?5?(20.4945?j1.1266)MVA 152?37.52?ST6?3??(0.386?j4.34)?(0.514?j5.7793)MVA35220.65052?0.54512?ST6?2??(0.386?j0.153)?(0.1345?j0.0533)MVA23526.3362?98.73692?ST6?1??(0.386?j6.584)?(3.2905?j56.1256)MVA352Sg?Sf??ST6?1?SG5??ST6?2??ST6?3?Si?35?(25.5114?j194.12)MVA計算每一個FLR的功率分布和電壓分布計算如下： FLR1：
?VT2?PR?QX40?1.346?45?31.7625??12.8970kVVN115 Vb?115??VT2?102.1030kVPR?QX10.0442?8.4?0.1442?16?VL1???0.8489kVVb102.1030Va?Vb??VL1?101.2541kVFLR2：
功率分布：
SL2?ZZ?Z*T3*L2*Sd?T3(Vb?VN)Z?ZL2****?VN?(0.092?j5.042)?(132.9792?j149.229)?1418.6727.392?j57.042T3?(4.8812?j13.8097)MVAST3?ZZ?Z*L2*L2*Sd?T3(Vb?VN)Z?ZL2?VN?(27.3?j52)?(132.9792?j149.229)?1418.6727.392?j57.042T3?(108.687?j122.62)MVA 電壓分布：
Sc1?SL2??SL2?(4.8812?j13.8097)?(2.67?j5.0854)?(7.5512?j8.7243)MVA7.5512?27.3?8.7243?52??2.424kV102.1030Vd?Vb??VL2?102.103?(?2.424)?104.527kV?VL2?功率分布：
FLR3：
SL3?ZZ?Z*T4*L3*Se?T4(VG3?Vd)Z?ZL3****?VN?(1.725?j39.703)?(89.945?j130.0151)?1037.9279.075?j53.73T4?(59.444?j16.846)MVAST4?ZZ?Z*L3*L3*Se?T4(Vb?VN)Z?ZL3?VN?(7.35?j14)?(89.945?j130.0151)?1037.9279.075?j53.73T4?(31.811?j60.1256)MVA 電壓分布：
Se1?SL3??SL3?(59.444?j19.846)?(24.51?j46.682)?(83.954?j26.836)MVA83.954?7.35?26.836?14?9.404kV105.5643Ve?Vd??VL3?96.16kV?VL3?功率分布：
FLR4：
SL4?ZZ?Z*T5*L4*Sf?T5(VG3?Vd)Z?ZL4****?VN=(0.368?j20.168)?(34.9449?j107.3469)?1037.92712.968?j44.168T5?(20.843?j19.689)MVAST4?ZZ?Z*L4*L4*Sf?T5(VG3?Vd)Z?ZL4?VN=(12.6?j24)?(34.9449?j107.3469)?1037.92712.968?j44.168T5?(1.398?j44.389)MVA 電壓分布：
Se1?SL3??SL3?(59.444?j16.846)?(24.51?j46.682)?(83.954?j63.528)MVA83.954?12.6?63.528?24?24.464kV105.5643Ve?Vd??VL3?81.10kV?VL4?FLR5: 這里我們先將f點和發(fā)電機G5當(dāng)做電源，經(jīng)過ZT61和ZT63構(gòu)成兩端供電網(wǎng)絡(luò)以g點作為運算負(fù)荷進行計算。ST6?ST4(0.386?j4.134)?(20.2656?j70.9293)?(22.0938?37)?35?(3.900?j25.1175)MVA0.772?j10.718(0.386?j6.584)?(20.2656?j70.9293)?(22.0938?37)?35??(16.5061?j91.7905)MVA0.772?j10.718電壓分布：
ST631?ST63??ST63?(16.6421?j97.5698)MVA16.6421?0.386?97.5698?4.134?10.9186kV37Vg?37??VT63?26.0814V?VT63?20.2656?0.386?70.9293?(?0.153)??0.1162kV
26.0814Vi?Vg??VT62?26.1976?VT62?20.4945?2.1?1.1266?4?1.815kV26.1976Vh?Vi??VL5?24.3826?VL5?
4設(shè)計圖及源程序
4.1MATLAB仿真
相關(guān)的原始數(shù)據(jù)輸入格式如下：
1、B1是支路參數(shù)矩陣，第一列和第二列是節(jié)點編號。節(jié)點編號由小到大編寫。
2、對于含有變壓器的支路，第一列為低壓側(cè)節(jié)點編號，第二列為高壓側(cè)節(jié)點編號，將變壓器的串聯(lián)阻抗置于低壓側(cè)處理，第三列為支路的串列阻抗參數(shù)，第四列為支路的對地導(dǎo)納參數(shù)，第五烈為含變壓器支路的變壓器的變比，第六列為變壓器是否是否含有變壓器的參數(shù)，其中“1”為含有變壓器，“0”為不含有變壓器。
3、B2為節(jié)點參數(shù)矩陣，其中第一列為節(jié)點注入發(fā)電功率參數(shù)；第二列為節(jié)點負(fù)荷功率參數(shù)；第三列為節(jié)點電壓參數(shù)；第六列為節(jié)點類型參數(shù)，其中“1”為平衡節(jié)點，“2”為PQ節(jié)點，“3”為PV節(jié)點參數(shù)。
4、X為節(jié)點號和對地參數(shù)矩陣。其中第一列為節(jié)點編號，第二列為節(jié)點對地參數(shù)。
4.2潮流計算源程序
%本程序的功能是用牛頓——拉夫遜法進行11節(jié)點潮流計算 clear;n=11;%input('請輸入節(jié)點數(shù)：n=');nl=11;%input('請輸入支路數(shù)：nl=');isb=1;%input('請輸入平衡母線節(jié)點號：isb=');pr=0.00001;%input('請輸入誤差精度：pr=');B1=[1
0.03512+0.08306i
0.13455i
0;
0.0068+0.18375i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.00811+0.24549i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.04215+0.09967i
0.04037i
0;
0.0068+0.18375i
0
1.02381
1;
0.02810+0.06645i
0.10764i
0;
0.05620+0.13289i
0.05382i
0;0.00811+0.24549i
0
1;
0.03512+0.08306i
0.13455i
0] B2=[0
0
1.1
1.1
0
1;
0
0
0
0
2;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0.204+0.12638i
0
0
2;
0
0
0
0
2;
0
0.306+0.18962i
0
0
2;
0
0
0
0
2;
0.5
0
1.1
1.1
0
3;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2];% B1矩陣：
1、支路首端號；
2、末端號；
3、支路阻抗；
4、支路對地電納 %
5、支路的變比；
6、支路首端處于K側(cè)為1，1側(cè)為0 % B2矩陣：
1、該節(jié)點發(fā)電機功率；
2、該節(jié)點負(fù)荷功率；
3、節(jié)點電壓初始值 %
4、PV節(jié)點電壓V的給定值；
5、節(jié)點所接的無功補償設(shè)備的容量 %
6、節(jié)點分類標(biāo)號：1為平衡節(jié)點（應(yīng)為1號節(jié)點）；2為PQ節(jié)點； %
3為PV節(jié)點；
%input('請輸入各節(jié)點參數(shù)形成的矩陣： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl
%支路數(shù)
if B1(i,6)==0
%左節(jié)點處于1側(cè)
p=B1(i,1);q=B1(i,2);
else
%左節(jié)點處于K側(cè)
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非對角元
Y(q,p)=Y(p,q);
%非對角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%對角元K側(cè)
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
%對角元1側(cè)
end %求導(dǎo)納矩陣
disp('導(dǎo)納矩陣 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);
%分解出導(dǎo)納陣的實部和虛部
for i=1:n
%給定各節(jié)點初始電壓的實部和虛部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
%PV節(jié)點電壓給定模值
end for i=1:n
%給定各節(jié)點注入功率
S(i)=B2(i,1)-B2(i,2);
%i節(jié)點注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5);%i節(jié)點無功補償量
end %===================== P=real(S);Q=imag(S);
%分解出各節(jié)點注入的有功和無功功率 ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次數(shù)ICT1、a；不滿足收斂要求的節(jié)點數(shù)IT2 while IT2~=0
% N0=2*n 雅可比矩陣的階數(shù)；N=N0+1擴展列
IT2=0;a=a+1;
for i=1:n
if i~=isb
%非平衡節(jié)點
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
P1=C(i)*e(i)+f(i)*D(i);%節(jié)點功率P計算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%節(jié)點功率Q計算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i節(jié)點有功和無功功率P',Q'的計算值
V2=e(i)^2+f(i)^2;
%電壓模平方
%========= 以下針對非PV節(jié)點來求取功率差及Jacobi矩陣元素 =========
if B2(i,6)~=3
%非PV節(jié)點
DP=P(i)-P1;
%節(jié)點有功功率差
DQ=Q(i)-Q1;
%節(jié)點無功功率差
%=============== 以上為除平衡節(jié)點外其它節(jié)點的功率計算 ================= %================= 求取Jacobi矩陣 ===================
for j1=1:n
if j1~=isb&j1~=i
%非平衡節(jié)點&非對角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de
X3=X2;
% X2=dp/df X3=dQ/de
X4=-X1;
% X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;J(p,N)=DQ;m=p+1;
% X3=dQ/de J(p,N)=DQ節(jié)點無功功率差
J(m,q)=X1;J(m,N)=DP;q=q+1;
% X1=dP/de J(m,N)=DP節(jié)點有功功率差
J(p,q)=X4;J(m,q)=X2;
% X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡節(jié)點&對角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%擴展列△Q
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%擴展列△P
J(m,q)=X2;
end
end
else
%=============== 下面是針對PV節(jié)點來求取Jacobi矩陣的元素 ===========
DP=P(i)-P1;
% PV節(jié)點有功誤差
DV=V(i)^2-V2;
% PV節(jié)點電壓誤差
for j1=1:n
if j1~=isb&j1~=i
%非平衡節(jié)點&非對角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
% dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i);
% dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV節(jié)點電壓誤差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV節(jié)點有功誤差
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡節(jié)點&對角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV節(jié)點電壓誤差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV節(jié)點有功誤差
J(m,q)=X2;
end
end
end
end
end %========= 以上為求雅可比矩陣的各個元素及擴展列的功率差或電壓差 =====================
for k=3:N0
% N0=2*n（從第三行開始，第一、二行是平衡節(jié)點）
k1=k+1;N1=N;
% N=N0+1 即 N=2*n+1擴展列△P、△Q 或 △U
for k2=k1:N1
% 從k+1列的Jacobi元素到擴展列的△P、△Q 或 △U
J(k,k2)=J(k,k2)./(J(k,k)+eps);% 用K行K列對角元素去除K行K列后的非對角元素進行規(guī)格化
end
J(k,k)=1;
% 對角元規(guī)格化K行K列對角元素賦1
%==================== 回代運算
=======================================
if k~=3
% 不是第三行
k > 3
k4=k-1;
for k3=3:k4
% 用k3行從第三行開始到當(dāng)前行的前一行k4行消去
for k2=k1:N1 %
k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當(dāng)前行k列元素消為0）
end
%用當(dāng)前行K2列元素減去當(dāng)前行k列元素乘以第k行K14 列元素
J(k3,k)=0;%當(dāng)前行第k列元素已消為0
end
if k==N0
%若已到最后一行
break;
end
%================== 前代運算
==================================
for k3=k1:N0
% 從k+1行到2*n最后一行
for k2=k1:N1
% 從k+1列到擴展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算
end
%用當(dāng)前行K2列元素減去當(dāng)前行k列元素乘以第k行K2列元素
J(k3,k)=0;%當(dāng)前行第k列元素已消為0
end
else
%是第三行k=3
%====================== 第三行k=3的前代運算 ========================
for k3=k1:N0
%從第四行到2n行（最后一行）
for k2=k1:N1
%從第四列到2n+1列（即擴展列）
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當(dāng)前行3列元素消為0）
end
%用當(dāng)前行K2列元素減去當(dāng)前行3列元素乘以第三行K2列元素
J(k3,k)=0;%當(dāng)前行第3列元素已消為0
end
end
end %====上面是用線性變換方式高斯消去法將Jacobi矩陣化成單位矩陣=====
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改節(jié)點電壓實部
k1=k+1;
f(L)=f(L)-J(k1,N);
%修改節(jié)點電壓虛部
end
%------修改節(jié)點電壓-----------
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
%電壓偏差量是否滿足要求
IT2=IT2+1;%不滿足要求的節(jié)點數(shù)加1
end
end
ICT2(a)=IT2;
%不滿足要求的節(jié)點數(shù)
ICT1=ICT1+1;
%迭代次數(shù) end %用高斯消去法解”w=-J*V“
disp('迭代次數(shù)：');disp(ICT1);disp('沒有達(dá)到精度要求的個數(shù)：');disp(ICT2);for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
%計算各節(jié)點電壓的模值
sida(k)=atan(f(k)./e(k))*180./pi;
%計算各節(jié)點電壓的角度
E(k)=e(k)+f(k)*j;
%將各節(jié)點電壓用復(fù)數(shù)表示 end %=============== 計算各輸出量 =========================== disp('各節(jié)點的實際電壓標(biāo)幺值E為(節(jié)點號從小到大排列)：');disp(E);
%顯示各節(jié)點的實際電壓標(biāo)幺值E用復(fù)數(shù)表示 disp('----------------------');disp('各節(jié)點的電壓大小V為(節(jié)點號從小到大排列)：');disp(V);
%顯示各節(jié)點的電壓大小V的模值 disp('----------------------');disp('各節(jié)點的電壓相角sida為(節(jié)點號從小到大排列)：');disp(sida);
%顯示各節(jié)點的電壓相角 for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q));%計算各節(jié)點的注入電流的共軛值
end
S(p)=E(p)*C(p);
%計算各節(jié)點的功率 S = 電壓 X 注入電流的共軛值 end disp('各節(jié)點的功率S為(節(jié)點號從小到大排列)：');disp(S);
%顯示各節(jié)點的注入功率
disp('----------------------');disp('各條支路的首端功率Si為(順序同您輸入B1時一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
end
disp(Si(p,q));
SSi(p,q)=Si(p,q);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
disp(ZF);
disp('----------------------');end disp('各條支路的末端功率Sj為(順序同您輸入B1時一致)：');
for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
end
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
disp(ZF);
disp('----------------------');end disp('各條支路的功率損耗DS為(順序同您輸入B1時一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
disp(ZF);
disp('----------------------');end
%本程序的功能是用牛頓——拉夫遜法進行10節(jié)點潮流計算 %本程序的功能是用牛頓——拉夫遜法進行潮流計算 clear;n=10;%input('請輸入節(jié)點數(shù)：n=');nl=10;%input('請輸入支路數(shù)：nl=');isb=1;%input('請輸入平衡母線節(jié)點號：isb=');pr=0.00001;%input('請輸入誤差精度：pr=');B1=[1
0.03512+0.08306i
0.13455i
0;
0.0068+0.18375i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.00811+0.24549i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.04215+0.09967i
0.04037i
0;
0.0068+0.18375i
0
1.02381
1;
0.02810+0.06645i
0.10764i
0;0.00811+0.24549i
0
1;
0.03512+0.08306i
0.13455i
0] B2=[0
0
1.1
1.1
0
1;
0
0
0
0
2;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0.204+0.12638i
0
0
2;
0
0
0
0
2;
0
0.306+0.18962i
0
0
2;
0
0
0
0
2;
0.5
0
1.1
1.1
0
3;
0
0.343+0.21256i
0
0
2];% B1矩陣：
1、支路首端號；
2、末端號；
3、支路阻抗；
4、支路對地電納 %
5、支路的變比；
6、支路首端處于K側(cè)為1，1側(cè)為0 % B2矩陣：
1、該節(jié)點發(fā)電機功率；
2、該節(jié)點負(fù)荷功率；
3、節(jié)點電壓初始值 %
4、PV節(jié)點電壓V的給定值；
5、節(jié)點所接的無功補償設(shè)備的容量 %
6、節(jié)點分類標(biāo)號：1為平衡節(jié)點（應(yīng)為1號節(jié)點）；2為PQ節(jié)點； %
3為PV節(jié)點；
%input('請輸入各節(jié)點參數(shù)形成的矩陣： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl
%支路數(shù)
if B1(i,6)==0
%左節(jié)點處于1側(cè)
p=B1(i,1);q=B1(i,2);
else
%左節(jié)點處于K側(cè)
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非對角元
Y(q,p)=Y(p,q);
%非對角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%對角元K側(cè)
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
%對角元1側(cè)
end %求導(dǎo)納矩陣
disp('導(dǎo)納矩陣 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);
%分解出導(dǎo)納陣的實部和虛部
for i=1:n
%給定各節(jié)點初始電壓的實部和虛部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
%PV節(jié)點電壓給定模值
end for i=1:n
%給定各節(jié)點注入功率
S(i)=B2(i,1)-B2(i,2);
%i節(jié)點注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5);%i節(jié)點無功補償量
end %===================== P=real(S);Q=imag(S);
%分解出各節(jié)點注入的有功和無功功率
ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次數(shù)ICT1、a；不滿足收斂要求的節(jié)點數(shù)IT2 while IT2~=0
% N0=2*n 雅可比矩陣的階數(shù)；N=N0+1擴展列
IT2=0;a=a+1;
for i=1:n
if i~=isb
%非平衡節(jié)點
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
P1=C(i)*e(i)+f(i)*D(i);%節(jié)點功率P計算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%節(jié)點功率Q計算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i節(jié)點有功和無功功率P',Q'的計算值
V2=e(i)^2+f(i)^2;
%電壓模平方
%========= 以下針對非PV節(jié)點來求取功率差及Jacobi矩陣元素 =========
if B2(i,6)~=3
%非PV節(jié)點
DP=P(i)-P1;
%節(jié)點有功功率差
DQ=Q(i)-Q1;
%節(jié)點無功功率差
%=============== 以上為除平衡節(jié)點外其它節(jié)點的功率計算 ================= %================= 求取Jacobi矩陣 ===================
for j1=1:n
if j1~=isb&j1~=i
%非平衡節(jié)點&非對角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de
X3=X2;
% X2=dp/df X3=dQ/de
X4=-X1;
% X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;J(p,N)=DQ;m=p+1;
% X3=dQ/de J(p,N)=DQ節(jié)點無功功率差
J(m,q)=X1;J(m,N)=DP;q=q+1;
% X1=dP/de J(m,N)=DP節(jié)點有功功率差
J(p,q)=X4;J(m,q)=X2;
% X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡節(jié)點&對角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%擴展列△Q
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%擴展列△P
J(m,q)=X2;
end
end
else
%=============== 下面是針對PV節(jié)點來求取Jacobi矩陣的元素
===========
DP=P(i)-P1;
% PV節(jié)點有功誤差
DV=V(i)^2-V2;
% PV節(jié)點電壓誤差
for j1=1:n
if j1~=isb&j1~=i
%非平衡節(jié)點&非對角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
% dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i);
% dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV節(jié)點電壓誤差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV節(jié)點有功誤差
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡節(jié)點&對角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV節(jié)點電壓誤差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV節(jié)點有功誤差
J(m,q)=X2;
end
end
end
end
end %========= 以上為求雅可比矩陣的各個元素及擴展列的功率差或電壓差 =====================
for k=3:N0
% N0=2*n（從第三行開始，第一、二行是平衡節(jié)點）
k1=k+1;N1=N;
% N=N0+1 即 N=2*n+1擴展列△P、△Q 或 △U
for k2=k1:N1
% 從k+1列的Jacobi元素到擴展列的△P、△Q 或 △U
J(k,k2)=J(k,k2)./J(k,k);% 用K行K列對角元素去除K行K列后的非對角元素進行規(guī)格化
end
J(k,k)=1;
% 對角元規(guī)格化K行K列對角元素賦1
%==================== 回代運算
=======================================
if k~=3
% 不是第三行
k > 3
k4=k-1;
for k3=3:k4
% 用k3行從第三行開始到當(dāng)前行的前一行k4行消
去
for k2=k1:N1 %
k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當(dāng)前行k列元素消為0）
end
%用當(dāng)前行K2列元素減去當(dāng)前行k列元素乘以第k行K2列元素
J(k3,k)=0;%當(dāng)前行第k列元素已消為0
end
if k==N0
%若已到最后一行
break;
end
%================== 前代運算
==================================
for k3=k1:N0
% 從k+1行到2*n最后一行
for k2=k1:N1
% 從k+1列到擴展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算
end
%用當(dāng)前行K2列元素減去當(dāng)前行k列元素乘以第k行K2列元素
J(k3,k)=0;%當(dāng)前行第k列元素已消為0
end
else
%是第三行k=3
%====================== 第三行k=3的前代運算 ========================
for k3=k1:N0
%從第四行到2n行（最后一行）
for k2=k1:N1
%從第四列到2n+1列（即擴展列）
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當(dāng)前行3列元素消為0）
end
%用當(dāng)前行K2列元素減去當(dāng)前行3列元素乘以第三行K2列元素
J(k3,k)=0;%當(dāng)前行第3列元素已消為0
end
end
end %====上面是用線性變換方式高斯消去法將Jacobi矩陣化成單位矩陣=====
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改節(jié)點電壓實部
k1=k+1;
f(L)=f(L)-J(k1,N);
%修改節(jié)點電壓虛部
end
%------修改節(jié)點電壓-----------
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
%電壓偏差量是否滿足要求
IT2=IT2+1;%不滿足要求的節(jié)點數(shù)加1
end
end
ICT2(a)=IT2;
%不滿足要求的節(jié)點數(shù)
ICT1=ICT1+1;
%迭代次數(shù) end %用高斯消去法解”w=-J*V“ disp('迭代次數(shù)：');disp(ICT1);disp('沒有達(dá)到精度要求的個數(shù)：');disp(ICT2);for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
%計算各節(jié)點電壓的模值
sida(k)=atan(f(k)./e(k))*180./pi;
%計算各節(jié)點電壓的角度
E(k)=e(k)+f(k)*j;
%將各節(jié)點電壓用復(fù)數(shù)表示 end %=============== 計算各輸出量 =========================== disp('各節(jié)點的實際電壓標(biāo)幺值E為(節(jié)點號從小到大排列)：');disp(E);
%顯示各節(jié)點的實際電壓標(biāo)幺值E用復(fù)數(shù)表示 disp('----------------------');disp('各節(jié)點的電壓大小V為(節(jié)點號從小到大排列)：');disp(V);
%顯示各節(jié)點的電壓大小V的模值 disp('----------------------');disp('各節(jié)點的電壓相角sida為(節(jié)點號從小到大排列)：');disp(sida);
%顯示各節(jié)點的電壓相角 for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q));%計算各節(jié)點的注入電流的共軛值
end
S(p)=E(p)*C(p);
%計算各節(jié)點的功率 S = 電壓 X 注入電流的共軛值 end disp('各節(jié)點的功率S為(節(jié)點號從小到大排列)：');disp(S);
%顯示各節(jié)點的注入功率
disp('----------------------');disp('各條支路的首端功率Si為(順序同您輸入B1時一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
end
disp(Si(p,q));
SSi(p,q)=Si(p,q);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
disp(ZF);
disp('----------------------');end disp('各條支路的末端功率Sj為(順序同您輸入B1時一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
end
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
disp(ZF);
disp('----------------------');end disp('各條支路的功率損耗DS為(順序同您輸入B1時一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
disp(ZF);
disp('----------------------');end
%本程序的功能是用牛頓——拉夫遜法進行12節(jié)點潮流計算 %本程序的功能是用牛頓——拉夫遜法進行潮流計算 clear;n=12;%input('請輸入節(jié)點數(shù)：n=');nl=12;%input('請輸入支路數(shù)：nl=');isb=1;%input('請輸入平衡母線節(jié)點號：isb=');pr=0.00001;%input('請輸入誤差精度：pr=');B1=[1
0.03512+0.08306i
0.13455i
0;
0.0068+0.18375i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.00811+0.24549i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.04215+0.09967i
0.04037i
0;
0.0068+0.18375i
0
1.02381
1;
0.02810+0.06645i
0.10764i
0;
0.05620+0.13289i
0.05382i
0;0.00811+0.24549i
0
1;
0.03512+0.08306i
0.13455i
0;
0.03512+0.08306i
0.13455i
0] B2=[0
0
1.1
1.1
0
1;
0
0
0
0
2;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0.204+0.12638i
0
0
2;
0
0
0
0
2;
0
0.306+0.18962i
0
0
2;
0
0
0
0
2;
0.5
0
1.1
1.1
0
3;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0
0
0
2];% B1矩陣：
1、支路首端號；
2、末端號；
3、支路阻抗；
4、支路對地電納 %
5、支路的變比；
6、支路首端處于K側(cè)為1，1側(cè)為0 % B2矩陣：
1、該節(jié)點發(fā)電機功率；
2、該節(jié)點負(fù)荷功率；
3、節(jié)點電壓初始值 %
4、PV節(jié)點電壓V的給定值；
5、節(jié)點所接的無功補償設(shè)備的容量 %
6、節(jié)點分類標(biāo)號：1為平衡節(jié)點（應(yīng)為1號節(jié)點）；2為PQ節(jié)點； %
3為PV節(jié)點；
%input('請輸入各節(jié)點參數(shù)形成的矩陣： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl
%支路數(shù)
if B1(i,6)==0
%左節(jié)點處于1側(cè)
p=B1(i,1);q=B1(i,2);
else
%左節(jié)點處于K側(cè)
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非對角元
Y(q,p)=Y(p,q);
%非對角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%對角元K側(cè)
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
%對角元1側(cè)
end %求導(dǎo)納矩陣
disp('導(dǎo)納矩陣 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);
%分解出導(dǎo)納陣的實部和虛部
for i=1:n
%給定各節(jié)點初始電壓的實部和虛部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
%PV節(jié)點電壓給定模值
end for i=1:n
%給定各節(jié)點注入功率
S(i)=B2(i,1)-B2(i,2);
%i節(jié)點注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5);%i節(jié)點無功補償量
end %===================== P=real(S);Q=imag(S);
%分解出各節(jié)點注入的有功和無功功率 ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次數(shù)ICT1、a；不滿足收斂要求的節(jié)點數(shù)IT2 while IT2~=0
% N0=2*n 雅可比矩陣的階數(shù)；N=N0+1擴展列
IT2=0;a=a+1;
for i=1:n
if i~=isb
%非平衡節(jié)點
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
P1=C(i)*e(i)+f(i)*D(i);%節(jié)點功率P計算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%節(jié)點功率Q計算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i節(jié)點有功和無功功率P',Q'的計算值
V2=e(i)^2+f(i)^2;
%電壓模平方
%========= 以下針對非PV節(jié)點來求取功率差及Jacobi矩陣元素 =========
if B2(i,6)~=3
%非PV節(jié)點
DP=P(i)-P1;
%節(jié)點有功功率差
DQ=Q(i)-Q1;
%節(jié)點無功功率差
%=============== 以上為除平衡節(jié)點外其它節(jié)點的功率計算 ================= %================= 求取Jacobi矩陣 ===================
for j1=1:n
if j1~=isb&j1~=i
%非平衡節(jié)點&非對角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de
X3=X2;
% X2=dp/df X3=dQ/de
X4=-X1;
% X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;J(p,N)=DQ;m=p+1;
% X3=dQ/de J(p,N)=DQ節(jié)點無功功率差
J(m,q)=X1;J(m,N)=DP;q=q+1;
% X1=dP/de J(m,N)=DP節(jié)點有功功率差
J(p,q)=X4;J(m,q)=X2;
% X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡節(jié)點&對角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%擴展列△Q
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%擴展列△P
J(m,q)=X2;
end
end
else
%=============== 下面是針對PV節(jié)點來求取Jacobi矩陣的元素 ===========
DP=P(i)-P1;
% PV節(jié)點有功誤差
DV=V(i)^2-V2;
% PV節(jié)點電壓誤差
for j1=1:n
if j1~=isb&j1~=i
%非平衡節(jié)點&非對角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
% dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i);
% dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV節(jié)點電壓誤差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV節(jié)點有功誤差
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡節(jié)點&對角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV節(jié)點電壓誤差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV節(jié)點有功誤差
J(m,q)=X2;
end
end
end
end
end %========= 以上為求雅可比矩陣的各個元素及擴展列的功率差或電壓差 =====================
for k=3:N0
% N0=2*n（從第三行開始，第一、二行是平衡節(jié)點）
k1=k+1;N1=N;
% N=N0+1 即 N=2*n+1擴展列△P、△Q 或 △U
for k2=k1:N1
% 從k+1列的Jacobi元素到擴展列的△P、△Q
或 △U
J(k,k2)=J(k,k2)./(J(k,k)+eps);% 用K行K列對角元素去除K行K列后的非對角元素進行規(guī)格化
end
J(k,k)=1;
% 對角元規(guī)格化K行K列對角元素賦1
%==================== 回代運算
=======================================
if k~=3
% 不是第三行
k > 3
k4=k-1;
for k3=3:k4
% 用k3行從第三行開始到當(dāng)前行的前一行k4行消去
for k2=k1:N1 %
k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當(dāng)前行k列元素消為0）
end
%用當(dāng)前行K2列元素減去當(dāng)前行k列元素乘以第k行K2列元素
J(k3,k)=0;%當(dāng)前行第k列元素已消為0
end
if k==N0
%若已到最后一行
break;
end
%================== 前代運算
==================================
for k3=k1:N0
% 從k+1行到2*n最后一行
for k2=k1:N1
% 從k+1列到擴展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算
end
%用當(dāng)前行K2列元素減去當(dāng)前行k列元素乘以第k行K2列元素
J(k3,k)=0;%當(dāng)前行第k列元素已消為0
end
else
%是第三行k=3
%====================== 第三行k=3的前代運算 ========================
for k3=k1:N0
%從第四行到2n行（最后一行）
for k2=k1:N1
%從第四列到2n+1列（即擴展列）
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當(dāng)前行3列元素消為0）
end
%用當(dāng)前行K2列元素減去當(dāng)前行3列元素乘以第三行K2列元素
J(k3,k)=0;%當(dāng)前行第3列元素已消為0
end
end
end %====上面是用線性變換方式高斯消去法將Jacobi矩陣化成單位矩陣=====
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改節(jié)點電壓實部
k1=k+1;
f(L)=f(L)-J(k1,N);
%修改節(jié)點電壓虛部
end
%------修改節(jié)點電壓-----------
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
%電壓偏差量是否滿足要求
IT2=IT2+1;%不滿足要求的節(jié)點數(shù)加1
end
end
ICT2(a)=IT2;
%不滿足要求的節(jié)點數(shù)
ICT1=ICT1+1;
%迭代次數(shù) end %用高斯消去法解”w=-J*V" disp('迭代次數(shù)：');disp(ICT1);disp('沒有達(dá)到精度要求的個數(shù)：');disp(ICT2);for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
%計算各節(jié)點電壓的模值
sida(k)=atan(f(k)./e(k))*180./pi;
%計算各節(jié)點電壓的角度
E(k)=e(k)+f(k)*j;
%將各節(jié)點電壓用復(fù)數(shù)表示 end %=============== 計算各輸出量 =========================== disp('各節(jié)點的實際電壓標(biāo)幺值E為(節(jié)點號從小到大排列)：');disp(E);
%顯示各節(jié)點的實際電壓標(biāo)幺值E用復(fù)數(shù)表示 disp('----------------------');disp('各節(jié)點的電壓大小V為(節(jié)點號從小到大排列)：');disp(V);
%顯示各節(jié)點的電壓大小V的模值 disp('----------------------');disp('各節(jié)點的電壓相角sida為(節(jié)點號從小到大排列)：');disp(sida);
%顯示各節(jié)點的電壓相角 for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q));%計算各節(jié)點的注入電流的共軛值
end
S(p)=E(p)*C(p);
%計算各節(jié)點的功率 S = 電壓 X 注入電流的共軛值 end disp('各節(jié)點的功率S為(節(jié)點號從小到大排列)：');disp(S);
%顯示各節(jié)點的注入功率
disp('----------------------');disp('各條支路的首端功率Si為(順序同您輸入B1時一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
end
disp(Si(p,q));
SSi(p,q)=Si(p,q);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
disp(ZF);
disp('----------------------');end disp('各條支路的末端功率Sj為(順序同您輸入B1時一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
end
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
disp(ZF);
disp('----------------------');end disp('各條支路的功率損耗DS為(順序同您輸入B1時一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
disp(ZF);
disp('----------------------');end
如果源程序的運行結(jié)果需要作圖可用下面的程序
figure(1);subplot(1,2,1);plot(V);xlabel('節(jié)點號');ylabel('電壓標(biāo)幺值');grid on;subplot(1,2,2);plot(sida);xlabel('節(jié)點號');ylabel('電壓角度');grid on;figure(2);subplot(2,2,1);P=real(S);Q=imag(S);bar(P);xlabel('節(jié)點號');ylabel('節(jié)點注入有功');grid on;subplot(2,2,2);bar(Q);xlabel('節(jié)點號');ylabel('節(jié)點注入無功');grid on;subplot(2,2,3);P1=real(Siz);Q1=imag(Siz);bar(P1);xlabel('支路號');ylabel('支路首端注入有功');grid on;subplot(2,2,4);bar(Q1);xlabel('支路號');ylabel('支路首端注入無功');grid on;
總結(jié)
通過本次課程設(shè)計讓我有復(fù)習(xí)了一次潮流計算的相關(guān)知識，跟家清晰了什么事潮流計算以及潮流計算的在電力系統(tǒng)的重要性。電力系統(tǒng)的穩(wěn)定運行狀況即是正常運行狀況，是指電力系統(tǒng)在穩(wěn)定運行條件下電壓、功率的分布，也稱為潮流分布。電力系統(tǒng)分析的潮流計算是電力系統(tǒng)分析的一個重要的部分。通過對電力系統(tǒng)潮流分布的分析和計算，可進一步對系統(tǒng)運行的安全性，經(jīng)濟性進行分析、評估，提出改進措施。同時潮流分布也是電力系統(tǒng)規(guī)劃設(shè)計的一項基礎(chǔ)工作。
整個計算過程的模型建立并不是十分復(fù)雜，但計算過程十分繁瑣、計算量相當(dāng)?shù)拇螅矣捎谥?jié)太多很容易算錯。不過在計算潮流計算的過程中卻對以往學(xué)過的電力系統(tǒng)分析的相關(guān)知識進行了一次較為深入的復(fù)習(xí)。而且整個計算對計算量的要求很大，鍛煉了我們的計算能力。而且對細(xì)節(jié)的把握也得到了鍛煉，做題的精細(xì)程度得到了提高。
參考文獻
[1]何仰贊, 溫增銀《電力系統(tǒng)分析》（第三版）［M].華中科技大學(xué)，2002 [2]http://baike.baidu.com/view/627420.[3]王守相，劉玉田 電力系統(tǒng)潮流計算研究現(xiàn)狀--《山東電力技術(shù)》1996年05期
[4]何仰贊，溫增銀.電力系統(tǒng)分析（上冊）第三版［M].湖北：華中科技大學(xué)出版社，2002 [5] 劉同娟.MATLAB在電路分析中的應(yīng)用.電氣電子教學(xué)學(xué)報.2002 [6] 西安交通大學(xué)等.電力系統(tǒng)計算[M].北京：水利電力出版社，1993.12 [7] 李光琦.電力系統(tǒng)暫態(tài)分析[M].北京： 水利電力出版社，2002.5 [8]何仰贊，溫增銀.電力系統(tǒng)分析（下冊）第三版［M].湖北：華中科技大學(xué)出版社，2002 [9]韋化，李濱，杭乃善，等．大規(guī)模水一火電力系統(tǒng)最優(yōu)潮流的現(xiàn)代內(nèi)點算法實現(xiàn)[J]．中國電機工程學(xué)報，2003．23(6)：13一l8．
[10]Chen Luo—nan，Suzuki Hideki，Katou Kazuo．Mean fieldtheory for optimal power flow[J]．IEEE Transactions OilPower Systems，1997，12(4)：1481·1486
第四篇：電力系統(tǒng)潮流計算
南 京 理 工 大 學(xué)
《電力系統(tǒng)穩(wěn)態(tài)分析》
課程報告
姓名
XX
學(xué) 號： 5*** 自動化學(xué)院 電氣工程
基于牛頓-拉夫遜法的潮流計算例題編程報學(xué)院(系): 專
業(yè): 題
目: 任課教師 碩士導(dǎo)師 告
楊偉 XX
2015年6月10號
基于牛頓-拉夫遜法的潮流計算例題編程報告
摘要：電力系統(tǒng)潮流計算的目的在于：確定電力系統(tǒng)的運行方式、檢查系統(tǒng)中各元件是否過壓或者過載、為電力系統(tǒng)繼電保護的整定提供依據(jù)、為電力系統(tǒng)的穩(wěn)定計算提供初值、為電力系統(tǒng)規(guī)劃和經(jīng)濟運行提供分析的基礎(chǔ)。潮流計算的計算機算法包含高斯—賽德爾迭代法、牛頓-拉夫遜法和P—Q分解法等，其中牛拉法計算原理較簡單、計算過程也不復(fù)雜，而且由于人們引入泰勒級數(shù)和非線性代數(shù)方程等在算法里從而進一步提高了算法的收斂性和計算速度。同時基于MATLAB的計算機算法以雙精度類型進行數(shù)據(jù)的存儲和運算, 數(shù)據(jù)精確度高，能進行潮流計算中的各種矩陣運算，使得傳統(tǒng)潮流計算方法更加優(yōu)化。
一 研究內(nèi)容
通過一道例題來認(rèn)真分析牛頓-拉夫遜法的原理和方法（采用極坐標(biāo)形式的牛拉法），同時掌握潮流計算計算機算法的相關(guān)知識，能看懂并初步使用MATLAB軟件進行編程，培養(yǎng)自己電力系統(tǒng)潮流計算機算法編程能力。
例題如下：用牛頓-拉夫遜法計算下圖所示系統(tǒng)的潮流分布，其中系統(tǒng)中5為平衡節(jié)點，節(jié)點5電壓保持U=1.05為定值，其他四個節(jié)點分別為PQ節(jié)點，給定的注入功率如圖所示。計算精度要求各節(jié)點電壓修正量不大于10-6。
二 牛頓-拉夫遜法潮流計算 1 基本原理
牛頓法是取近似解x(k)之后，在這個基礎(chǔ)上，找到比x(k)更接近的方程的根，一步步地迭代，找到盡可能接近方程根的近似根。牛頓迭代法其最大優(yōu)點是在方程f(x)=0的單根附近時誤差將呈平方減少，而且該法還可以用來求方程的重根、復(fù)根。電力系統(tǒng)潮流計算，一般來說，各個母線所供負(fù)荷的功率是已知的，各個節(jié)點的電壓是未知的（平衡節(jié)點外）可以根據(jù)網(wǎng)絡(luò)結(jié)構(gòu)形成節(jié)點導(dǎo)納矩陣，然后由節(jié)點導(dǎo)納矩陣列寫功率方程，由于功率方程里功率是已知的，電壓的幅值和相角是未知的，這樣潮流計算的問題就轉(zhuǎn)化為求解非線性方程組的問題了。為了便于用迭代法解方程組，需要將上述功率方程改寫成功率平衡方程，并對功率平衡方程求偏導(dǎo)，得出對應(yīng)的雅可比矩陣，給未知節(jié)點賦電壓初值，將初值帶入功率平衡方程，得到功率不平衡量，這樣由功率不平衡量、雅可比矩陣、節(jié)點電壓不平衡量（未知的）構(gòu)成了誤差方程，解誤差方程，得到節(jié)點電壓不平衡量，節(jié)點電壓加上節(jié)點電壓不平衡量構(gòu)成節(jié)點電壓新的初值，將新的初值帶入原來的功率平衡方程，并重新形成雅可比矩陣，然后計算新的電壓不平衡量，這樣不斷迭代，不斷修正，一般迭代三到五次就能收斂。2 基本步驟和設(shè)計流程圖
形成了雅克比矩陣并建立了修正方程式，運用牛頓-拉夫遜法計算潮流的核心問題已經(jīng)解決，已有可能列出基本計算步驟并編制流程圖。由課本總結(jié)基本步驟如下：
1)形成節(jié)點導(dǎo)納矩陣Y；
2)設(shè)各節(jié)點電壓的初值，如果是直角坐標(biāo)的話設(shè)電壓的實部e和虛部f；如果是極坐標(biāo)的話則設(shè)電壓的幅值U和相角a；
3)將各個節(jié)點電壓的初值代入公式求修正方程中的不平衡量以及修正方程的系數(shù)矩陣的雅克比矩陣；
4)解修正方程式，求各節(jié)點電壓的變化量，即修正量； 5)計算各個節(jié)點電壓的新值，即修正后的值；
6)利用新值從第（3）步開始進入下一次迭代，直至達(dá)到精度退出循環(huán)； 7)計算平衡節(jié)點的功率和線路功率，輸出最后計算結(jié)果； ① 公式推導(dǎo)
② 流程圖
三
matlab編程代碼
clear;
% 如圖所示1，2，3，4為PQ節(jié)點,5為平衡節(jié)點
y=0;
% 輸入原始數(shù)據(jù)，求節(jié)點導(dǎo)納矩陣
y(1,2)=1/(0.07+0.21j);
y(4,5)=0;y(1,3)=1/(0.06+0.18j);
y(1,4)=1/(0.05+0.10j);
y(1,5)=1/(0.04+0.12j);
y(2,3)=1/(0.05+0.10j);
y(2,5)=1/(0.08+0.24j);
y(3,4)=1/(0.06+0.18j);
for i=1:5
for j=i:5
y(j,i)=y(i,j);
end
end
Y=0;
% 求節(jié)點導(dǎo)納矩陣中互導(dǎo)納
for i=1:5
for j=1:5
if i~=j
Y(i,j)=-y(i,j);
end
end
end
% 求節(jié)點導(dǎo)納矩陣中自導(dǎo)納
for i=1:5
Y(i,i)=sum(y(i,:));
end
Y
% Y為導(dǎo)納矩陣
G=real(Y);
B=imag(Y);% 輸入原始節(jié)點的給定注入功率
S(1)=0.3+0.3j;
S(2)=-0.5-0.15j;
S(3)=-0.6-0.25j;
S(4)=-0.7-0.2j;
S(5)=0;
P=real(S);
Q=imag(S);
% 賦初值，U為節(jié)點電壓的幅值，a為節(jié)點電壓的相位角
U=ones(1,5);
U(5)=1.05;
a=zeros(1,5);
x1=ones(8,1);
x2=ones(8,1);
k=0;
while max(x2)>1e-6
for i=1:4
for j=1:4
H(i,j)=0;
N(i,j)=0;
M(i,j)=0;
L(i,j)=0;
oP(i)=0;
oQ(i)=0;
end
end
% 求有功、無功功率不平衡量
for i=1:4
for j=1:5
oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));
end
oP(i)=oP(i)+P(i);
oQ(i)=oQ(i)+Q(i);
end
x2=[oP,oQ]';
% x2為不平衡量列向量
% 求雅克比矩陣
% 當(dāng)i~=j時，求H,N,M,L
for i=1:4
for j=1:4
if i~=j
H(i,j)=-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));
N(i,j)=-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
L(i,j)=H(i,j);
M(i,j)=-N(i,j);
end
end
end
% 當(dāng)i=j時，求H,N,M,L
for i=1:4
for j=1:5
if i~=j H(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));N(i,i)=N(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))
end
end
N(i,i)=N(i,i)-2*(U(i))^2*G(i,i);
L(i,i)=L(i,i)+2*(U(i))^2*B(i,i);
end
J=[H,N;M,L]
% J為雅克比矩陣
x1=-((inv(J))*x2);
% x1為所求△x的列向量
% 求節(jié)點電壓新值，準(zhǔn)備下一次迭代
for i=1:4
oa(i)=x1(i);
oU(i)=x1(i+4)*U(i);
end
for i=1:4
a(i)=a(i)+oa(i);
U(i)=U(i)+oU(i);
end
k=k+1;
end
k,U,a
% 求節(jié)點注入功率
i=5;
for j=1:5
P(i)=U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)))+P(i);
Q(i)=U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))+Q(i);
end
S(5)=P(5)+Q(5)*sqrt(-1);
S
% 求節(jié)點注入電流
I=Y*U'
四
運行結(jié)果
節(jié)點導(dǎo)納矩陣
經(jīng)過五次迭代后的雅克比矩陣
迭代次數(shù)以及節(jié)點電壓的幅值和相角（弧度數(shù)）
節(jié)點注入功率和電流
五 結(jié)果分析
在這次學(xué)習(xí)和實際操作過程里：首先，對電力系統(tǒng)分析中潮流計算的部分特別是潮流計算的計算機算法中的牛頓-拉夫遜法進行深入的研讀，弄明白了其原理、計算過程、公式推導(dǎo)以及設(shè)計流程。牛頓-拉夫遜法是求解非線性方程的迭代過程，其計算公式為?F?J?X，式中J為所求函數(shù)的雅可比矩陣；?X為需要求的修正值；?F為不平衡的列向量。利用x(*)=x(k+1)+?X(k+1)進行多次迭代，通過迭代判據(jù)得到所需要的精度值即準(zhǔn)確值x(*)。六 結(jié)論
通過這個任務(wù)，自己在matlab編程，潮流計算，word文檔的編輯功能等方面均有提高，但也暴漏出一些問題：理論知識儲備不足，對matlab的性能和特點還不能有一個全面的把握，對word軟件也不是很熟練，相信通過以后的學(xué)習(xí)能彌補這些不足，達(dá)到一個新的層次。
第五篇：基于MATLAB的電力系統(tǒng)潮流計算
基于MATLAB的電力系統(tǒng)潮流計算
%簡單潮流計算的小程序，相關(guān)的原始數(shù)據(jù)數(shù)據(jù)數(shù)據(jù)輸入格式如下：
%B1是支路參數(shù)矩陣，第一列和第二列是節(jié)點編號。節(jié)點編號由小到大編寫 %對于含有變壓器的支路，第一列為低壓側(cè)節(jié)點編號，第二列為高壓側(cè)節(jié)點 %編號，將變壓器的串聯(lián)阻抗置于低壓側(cè)處理。%第三列為支路的串列阻抗參數(shù)。%第四列為支路的對地導(dǎo)納參數(shù)。
%第五烈為含變壓器支路的變壓器的變比
%第六列為變壓器是否是否含有變壓器的參數(shù)，其中“1”為含有變壓器，%“0”為不含有變壓器。
%B2為節(jié)點參數(shù)矩陣，其中第一列為節(jié)點注入發(fā)電功率參數(shù)；第二列為節(jié)點 %負(fù)荷功率參數(shù)；第三列為節(jié)點電壓參數(shù)；第六列為節(jié)點類型參數(shù)，其中 %“1”為平衡節(jié)點，“2”為PQ節(jié)點，“3”為PV節(jié)點參數(shù)。
%X為節(jié)點號和對地參數(shù)矩陣。其中第一列為節(jié)點編號，第二列為節(jié)點對地 %參數(shù)。
n=input('請輸入節(jié)點數(shù):n=');n1=input('請輸入支路數(shù):n1=');isb=input('請輸入平衡節(jié)點號:isb=');pr=input('請輸入誤差精度:pr=');B1=input('請輸入支路參數(shù):B1=');B2=input('請輸入節(jié)點參數(shù):B2=');X=input('節(jié)點號和對地參數(shù):X=');Y=zeros(n);Times=1;%置迭代次數(shù)為初始值 %創(chuàng)建節(jié)點導(dǎo)納矩陣 for i=1:n1 if B1(i,6)==0 %不含變壓器的支路 p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else %含有變壓器的支路 p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3);Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));end end Y OrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);%將OrgS、DetaS初始化
%創(chuàng)建OrgS，用于存儲初始功率參數(shù) h=0;j=0;for i=1:n %對PQ節(jié)點的處理 if i~=isb&B2(i,6)==2 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end for i=1:n %對PV節(jié)點的處理，注意這時不可再將h初始化為0 if i~=isb&B2(i,6)==3 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end OrgS %創(chuàng)建PVU 用于存儲PV節(jié)點的初始電壓 PVU=zeros(n-h-1,1);t=0;for i=1:n if B2(i,6)==3 t=t+1;PVU(t,1)=B2(i,3);end end PVU %創(chuàng)建DetaS，用于存儲有功功率、無功功率和電壓幅值的不平衡量 h=0;for i=1:n %對PQ節(jié)點的處理 if i~=isb&B2(i,6)==2 h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);end end t=0;for i=1:n %對PV節(jié)點的處理，注意這時不可再將h初始化為0 if i~=isb&B2(i,6)==3 h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;end end DetaS %創(chuàng)建I，用于存儲節(jié)點電流參數(shù) i=zeros(n-1,1);h=0;for i=1:n if i~=isb h=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));end end I %創(chuàng)建Jacbi(雅可比矩陣)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %對PQ節(jié)點的處理 if B2(i,6)==2 h=h+1;for j=1:n if j~=isb k=k+1;if i==j %對角元素的處理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else %非對角元素的處理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);end if k==(n-1)%將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行
k=0;end end end end end k=0;for i=1:n %對PV節(jié)點的處理 if B2(i,6)==3 h=h+1;for j=1:n if j~=isb k=k+1;if i==j %對角元素的處理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else %非對角元素的處理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;end if k==(n-1)%將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行
k=0;end end end end end Jacbi %求解修正方程，獲取節(jié)點電壓的不平衡量 DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU %修正節(jié)點電壓 j=0;for i=1:n %對PQ節(jié)點處理 if B2(i,6)==2 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end for i=1:n %對PV節(jié)點的處理 if B2(i,6)==3 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end B2 %開始循環(huán)********************************************************************** while abs(max(DetaU))>pr OrgS=zeros(2*n-2,1);%！!初始功率參數(shù)在迭代過程中是不累加的，所以在這里必須將其初始化為零矩陣 h=0;j=0;for i=1:n if i~=isb&B2(i,6)==2 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end for i=1:n if i~=isb&B2(i,6)==3 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end OrgS %創(chuàng)建DetaS h=0;for i=1:n if i~=isb&B2(i,6)==2 h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);end end t=0;for i=1:n if i~=isb&B2(i,6)==3 h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;end end DetaS %創(chuàng)建I i=zeros(n-1,1);h=0;for i=1:n if i~=isb h=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));end end I %創(chuàng)建Jacbi Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n if B2(i,6)==2 h=h+1;for j=1:n if j~=isb k=k+1;if i==j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);end if k==(n-1)k=0;end end end end end k=0;for i=1:n if B2(i,6)==3 h=h+1;for j=1:n if j~=isb k=k+1;if i==j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;end if k==(n-1)k=0;end end end end end Jacbi DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU %修正節(jié)點電壓 j=0;for i=1:n if B2(i,6)==2 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end for i=1:n if B2(i,6)==3 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end B2 Times=Times+1;%迭代次數(shù)加1 end Times 一個原始數(shù)據(jù)的例子 節(jié)點數(shù) 5 支路數(shù) 5平衡節(jié)點編號 5 精度pr 0.000001 B1(支路參數(shù)矩陣)[1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 3 0.03i 0 1.05 1] B2(節(jié)點參數(shù)矩陣)[0-1.6-0.8i 1 0 0 2;0-2-1i 1 0 0 2;0-3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1] X(節(jié)點號和對地參數(shù))[1 0;2 0;3 0;4 0;5 0]
電力系統(tǒng)潮流計算
——9結(jié)點算例-PQ法
原始數(shù)據(jù)錄入data.txt文檔:
標(biāo)號,起始結(jié)點,終止結(jié)點,支路電阻參數(shù),支路電抗參數(shù),支路對地導(dǎo)納參數(shù) 1,2,5,0.0,0.063,0.0, 2,5,9,0.019,0.072,0.075, 3,6,9,0.012,0.101,0.105, 4,3,6,0.0,0.059,0.0, 5,6,8,0.039,0.17,0.179, 6,4,8,0.017,0.092,0.079, 7,5,7,0.032,0.161,0.153, 8,4,7,0.01,0.085,0.088, 9,1,4,0.0,0.058,0.0, 潮流程序chaoliu.txt文檔: #include ; float x1[N-1],x2 ;for(i=1;i guass(1,N-1,y1,B,x1);for(i=1;i guass(N-M,M,y2,B,x2);for(i=N-M;i else { kp=0;if(kq==0)val(u,g,b,r,ku,kr,h);else goto top;} } } void val(float u[N],float g[N][N],float b[N][N],float r[N],int ku, int kr,float h[N][N]){ float ps=0,pv1=0,pv2=0;float qs=0,qv1=0,qv2=0;float p[N][N]={0};float q[N][N]={0};float s[N][N];float dp[N][N]={0};float dq[N][N]={0};float ds[N][N];float dSp=0,dSq=0;int i,j;FILE *fp1;printf(”n=====ping heng jie dian gong lv =====n“);getch();for(i=0;i printf(”n=======shu ju bao cun=====n“);fp1=fopen(”jieguo.txt“,”w+“);{ fprintf(fp1,”xian lu cao liu:n“);for(i=0;i
點擊咨詢 