标签：分析 end 暂态 linedata Vm nl nr Ybus 电力系统
http://www.apollocode.net/a/1086.html
电力系统暂态稳定分析程序，
clear
clc
basemva = 100; accuracy = 0.0001; maxiter = 10;
% Bus BusbVoltage Angle ---Load---- -----Generator-----Static Mvar
% No code Mag. Degree MW Mvar MW Mvar Qmin Qmax Qc/-Ql
busdata=[1 1 1.06 0.0 00.00 00.00 0.00 00.00 0 0 0
2 2 1.04 0.0 00.00 00.00 150.00 00.00 0 140 0
3 2 1.03 0.0 00.00 00.00 100.00 00.00 0 90 0
4 0 1.0 0.0 100.00 70.00 00.00 00.00 0 0 0
5 0 1.0 0.0 90.00 30.00 00.00 00.00 0 0 0
6 0 1.0 0.0 160.00 110.00 00.00 00.00 0 0 0];
% Line code
% Bus bus R X 1/2 B = 1 for lines
% nl nr p.u. p.u. p.u. >1 or<1 tr. tap at bus nl
linedata=[1 4 0.035 0.225 0.0065 1.0
1 5 0.025 0.105 0.0045 1.0
1 6 0.040 0.215 0.0055 1.0
2 4 0.000 0.035 0.0000 1.0
3 5 0.000 0.042 0.0000 1.0
4 6 0.028 0.125 0.0035 1.0
5 6 0.026 0.175 0.0300 1.0];
% Gen. Ra Xd' H
gendata=[ 1 0 0.20 20
2 0 0.15 4
3 0 0.25 5];
%%%-------------------- Ene of input data --------------------------------
%%------------------------Start creat Ybus----------------------------------
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; %branch admittance
for n = 1:nbr
if a(n) <= 0 a(n) = 1;
else end
Ybus=zeros(nbus,nbus);
for k=1:nbr;
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
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
%%------------------------END creat Ybus----------------------------------
%%------------------------Start Newton-Raphson load Flow ----------------------------------
ns=0; ng=0; Vm=0; delta=0; yload=0; deltad=0;
nbus = length(busdata(:,1));
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;
% Start of iterations
clear A DC J DX
while maxerror >= accuracy & iter <= maxiter
for i=1:m
for k=1:m
A(i,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 i=1:nbr
if nl(i) == n | nr(i) == n
if nl(i) == n, l = nr(i); end
if nr(i) == n, l = nl(i); 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);
A(nn, ll) =-Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));
if kb(l) == 0
A(nn, lk) =Vm(n)*Ym(n,l)*cos(t(n,l)- delta(n) + delta(l));end
if kb(n) == 0
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
A(lm, lk) =-Vm(n)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));end
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
if kb(n) == 2 Q(n)=Qk;
if Qmax(n) ~= 0
Qgc = Q(n)*basemva + Qd(n) - Qsh(n);
if iter <= 7
if iter > 2
if Qgc < Qmin(n),
Vm(n) = Vm(n) + 0.01;
elseif Qgc > Qmax(n),
Vm(n) = Vm(n) - 0.01;end
else, end
else,end
else,end
end
if kb(n) ~= 1
A(nn,nn) = J11;
DC(nn) = P(n)-Pk;
end
if kb(n) == 0
A(nn,lm) = 2*Vm(n)*Ym(n,n)*cos(t(n,n))+J22;
A(lm,nn)= J33;
A(lm,lm) =-2*Vm(n)*Ym(n,n)*sin(t(n,n))-J44;
DC(lm) = Q(n)-Qk;
end
end
DX=A\DC';
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.\n\n')
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);
elseif kb(n) ==2
k=k+1;
S(n)=P(n)+j*Q(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);
%%------------------------End Newton-Raphson load Flow ----------------------------------
%%------------------------Start transient Stability program ----------------------------------
global Pm f H E Y th ngg
f=60;
ngr=gendata(:,1);
ngg=length(gendata(:,1));
%%
for k=1:ngg
zdd(ngr(k))=gendata(k, 2)+j*gendata(k,3);
H(k)=gendata(k,4);
end
%%
for k=1:ngg
I=conj(S(ngr(k)))/conj(V(ngr(k)));
Ep(k) = V(ngr(k))+zdd(ngr(k))*I;
Pm(k)=real(S(ngr(k)));
end
E=abs(Ep); d0=angle(Ep);
for k=1:ngg
nl(nbr+k) = nbus+k;
nr(nbr+k) = gendata(k, 1);
R(nbr+k) = real(zdd(ngr(k)));
X(nbr+k) = imag(zdd(ngr(k)));
Bc(nbr+k) = 0;
a(nbr+k) = 1.0;
yload(nbus+k)=0;
end
nbr1=nbr; nbus1=nbus;
nbrt=nbr+ngg;
nbust=nbus+ngg;
linedata=[nl, nr, R, X, -j*Bc, a];
%%%%%%%%---------- creat Ybus before Fault -----------------------
%%-----------------Lowd flow YBUS -------------------------------
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;
for n = 1:nbr
if a(n) <= 0 a(n) = 1; else end
Ybus=zeros(nbus,nbus);
for k=1:nbr;
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
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
%%-----------------End Lowd flow YBUS -------------------------------
for k=1:nbust
Ybus(k,k)=Ybus(k,k)+yload(k);
end
YLL=Ybus(1:nbus1, 1:nbus1);
YGG = Ybus(nbus1+1:nbust, nbus1+1:nbust);
YLG = Ybus(1:nbus1, nbus1+1:nbust);
Ybf=YGG-YLG.'*inv(YLL)*YLG;
clear YLL YGG YLG
%%%%%%%%---------- END Ybus before Fault -----------------------
Y=abs(Ybf); th=angle(Ybf);
Pm=zeros(1, ngg);
for ii = 1:ngg
for jj = 1:ngg
Pm(ii) = Pm(ii) + E(ii)*E(jj)*Y(ii, jj)*cos(th(ii, jj)-d0(ii)+d0(jj));
end,
end
nf=input('Please Enter fault location (Bus Number ?) ');
%%%------------- Ybus during fault ---------------------------------------
nbusf=nbust-1;
Ybus(:,nf:nbusf)=Ybus(:,nf+1:nbust);
Ybus(nf:nbusf,:)=Ybus(nf+1:nbust,:);
YLL=Ybus(1:nbus1-1, 1:nbus1-1);
YGG = Ybus(nbus1:nbusf, nbus1:nbusf);
YLG = Ybus(1:nbus1-1, nbus1:nbusf);
Ydf=YGG-YLG.'*inv(YLL)*YLG;
%%%%%--------------------------------------------------------------
%%%%%--------------- Ybus after Fault -------------------------------
nl=linedata(:, 1);
nr=linedata(:, 2);
remove = 0;
while remove ~= 1
fline=input('Enter removed line? for example [5,6]---->? ');
nlf=fline(1); nrf=fline(2);
for k=1:nbrt
if nl(k)==nlf & nr(k)==nrf
remove = 1;
m=k;
else, end
end
if remove ~= 1
fprintf('\nThe line to be removed does not exist in the line data. try again.\n\n')
end
end
linedat2(1:m-1,:)= linedata(1:m-1,:);
linedat2(m:nbrt-1,:)=linedata(m+1:nbrt,:);
linedat0=linedata;
linedata=linedat2;
%%%%%-----------------Load Flow YBUS -------------------------------------
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;
for n = 1:nbr
if a(n) <= 0 a(n) = 1; else end
Ybus=zeros(nbus,nbus);
for k=1:nbr;
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
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
%%%%%-----------------End of Load Flow YBUS -------------------------------------
for k=1:nbust
Ybus(k,k)=Ybus(k,k)+yload(k);
end
YLL=Ybus(1:nbus1, 1:nbus1);
YGG = Ybus(nbus1+1:nbust, nbus1+1:nbust);
YLG = Ybus(1:nbus1, nbus1+1:nbust);
Yaf=YGG-YLG.'*inv(YLL)*YLG;
linedata=linedat0;
%%%%%--------------- End Ybus after Fault -------------------------------
tc=input('Enter time of removing fault in sec. tc = ');
tf=input('Enter simulation time in sec. tf = ');
clear t x del
t0 = 0;
w0=zeros(1, length(d0));
x0 = [d0, w0];
tol=0.0001;
Y=abs(Ydf); th=angle(Ydf);
tspan=[t0, tc];
[t1, xf] =ode23('trstability', tspan, x0);
x0c =xf(length(xf), :);
Y=abs(Yaf); th=angle(Yaf);
tspan = [tc, tf];
[t2,xc] =ode23('trstability', tspan, x0c);
t =[t1; t2]; x = [xf; xc];
for k=1:nbus1
if kb(k)==1
ms=k;
end
end
kk=0;
for k=1:ngg
if k~=ms
kk=kk+1;
del(:,kk)=180/pi*(x(:,k)-x(:,ms));
else, end
end
h=figure; figure(h)
plot(t, del)
title(['Phase angle difference (fault cleared at ', num2str(tc),'s)'])
xlabel('t, sec'), ylabel('Delta, degree'), grid
%%------------------------End transient Stability program ----------------------------------
标签：分析,end,暂态,linedata,Vm,nl,nr,Ybus,电力系统 来源： https://blog.csdn.net/ucas_123/article/details/116106995
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。