coding new england 39-bus system

6 views (last 30 days)
Sulaiman AlEidan
Sulaiman AlEidan on 7 Apr 2015
Answered: rajesh yerupalli on 1 Jul 2019
first of all this is the code
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); % new end for k=1:ngg I=conj(S(ngr(k)))/conj(V(ngr(k))); Ep(k) = V(ngr(k))+zdd(ngr(k))*I; % new Pm(k)=real(S(ngr(k))); % new 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]; [Ybus, Ybf]=ybusbf(linedata, yload, nbus1,nbust); fprintf('\nPrefault reduced bus admittance matrix \n') Ybf Y=abs(Ybf); th=angle(Ybf); Pm=zeros(1, ngg); disp([' G(i) E''(i) d0(i) Pm(i)']) 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, fprintf(' %g', ngr(ii)), fprintf(' %8.4f',E(ii)), fprintf(' %8.4f', 180/pi*d0(ii)) fprintf(' %8.4f \n',Pm(ii)) end respfl='y'; while respfl =='y' respfl=='Y' nf=input('Enter faulted bus No. -> '); rtn=isempty(nf); if rtn==1; nf=-1; end while nf 0 nf > nbus fprintf('Faulted bus No. must be between 1 & %g \n', nbus) nf = input('Enter Faulted Bus No. - '); rtn=isempty(nf); if rtn==1; nf=-1; end end
fprintf('\nFaulted reduced bus admittance matrix\n') Ydf=ybusdf(Ybus, nbus1, nbust, nf) %Fault cleared [Yaf]=ybusaf(linedata, yload, nbus1,nbust, nbrt); fprintf('\nPostfault reduced bus admittance matrix\n') Yaf resptc='y'; while resptc =='y' resptc=='Y' tc=input('Enter clearing time of fault in sec. tc = '); tf=input('Enter final 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); %[t1, xf] =ode23('dfpek', t0, tc, x0, tol); % Solution during fault (use with MATLAB 4) tspan=[t0, tc]; %use with MATAB 5 [t1, xf] =ode23('dfpek', tspan, x0); % Solution during fault (use with MATLAB 5) x0c =xf(length(xf), :); Y=abs(Yaf); th=angle(Yaf); %[t2,xc] =ode23('afpek', tc, tf, x0c, tol); % Postfault solution (use with MATLAB 4) tspan = [tc, tf]; % use with MATLAB 5 [t2,xc] =ode23('afpek', tspan, x0c); % Postfault solution (use with MATLAB 5) t =[t1; t2]; x = [xf; xc]; fprintf('\nFault is cleared at %4.3f Sec. \n', tc) for k=1:nbus if kb(k)==1 ms=k; else, end end fprintf('\nPhase angle difference of each machine \n') fprintf('with respect to the slack in degree.\n') fprintf(' t - sec') kk=0; for k=1:ngg if k~=ms kk=kk+1; del(:,kk)=180/pi*(x(:,k)-x(:,ms)); fprintf(' d(%g,',ngr(k)), fprintf('%g)', ngr(ms)) else, end end fprintf(' \n') disp([t, del]) h=figure; figure(h) plot(t, del) title(['Phase angle difference (fault cleared at ', num2str(tc),'s)']) xlabel('t, sec'), ylabel('Delta, degree'), grid resp=0; while strcmp(resp, 'n')~=1 && strcmp(resp, 'N')~=1 && strcmp(resp, 'y')~=1 && strcmp(resp, 'Y')~=1 resp=input('Another clearing time of fault? Enter ''y'' or ''n'' within quotes -> ');
if strcmp(resp, 'n')~=1 && strcmp(resp, 'N')~=1 && strcmp(resp, 'y')~=1 && strcmp(resp, 'Y')~=1
fprintf('\n Incorrect reply, try again \n\n'), end
end
resptc=resp;
end resp2=0; while strcmp(resp2, 'n')~=1 && strcmp(resp2, 'N')~=1 && strcmp(resp2, 'y')~=1 && strcmp(resp2, 'Y')~=1 resp2=input('Another fault location: Enter ''y'' or ''n'' within quotes -> '); if strcmp(resp2, 'n')~=1 && strcmp(resp2, 'N')~=1 && strcmp(resp2, 'y')~=1 && strcmp(resp2, 'Y')~=1 fprintf('\n Incorrect reply, try again \n\n'), end respf1=resp2; end if respf1=='n' respf1=='N', return, else, end end
and I wanna simulate this system to have different clearing time and fault location
clear all basemva = 100; accuracy = 0.0001; maxiter = 100;
busdata=[1 0 1 0 0.0 0.0 0.0 0.0 0 0 0 2 0 1 0 0.0 0.0 0.0 0.0 0 0 0 3 0 1 0 322.0 2.4 0.0 0.0 0 0 0 4 0 1 0 500.0 184.0 0.0 0.0 0 0 0 5 0 1 0 0.0 0.0 0.0 0.0 0 0 0 6 0 1 0 0.0 0.0 0.0 0.0 0 0 0 7 0 1 0 233.8 84.0 0.0 0.0 0 0 0 8 0 1 0 522.0 176.0 0.0 0.0 0 0 0 9 0 1 0 0.0 0.0 0.0 0.0 0 0 0 10 0 1 0 0.0 0.0 0.0 0.0 0 0 0 11 0 1 0 0.0 0.0 0.0 0.0 0 0 0 12 0 1 0 7.5 88.0 0.0 0.0 0 0 0 13 0 1 0 0.0 0.0 0.0 0.0 0 0 0 14 0 1 0 0.0 0.0 0.0 0.0 0 0 0 15 0 1 0 320.0 153.0 0.0 0.0 0 0 0 16 0 1 0 329.0 32.3 0.0 0.0 0 0 0 17 0 1 0 0.0 0.0 0.0 0.0 0 0 0 18 0 1 0 158.0 30.0 0.0 0.0 0 0 0 19 0 1 0 0.0 0.0 0.0 0.0 0 0 0 20 0 1 0 628.0 103.0 0.0 0.0 0 0 0 21 0 1 0 274.0 115.0 0.0 0.0 0 0 0 22 0 1 0 0.0 0.0 0.0 0.0 0 0 0 23 0 1 0 247.5 84.6 0.0 0.0 0 0 0 24 0 1 0 308.6 -92.0 0.0 0.0 0 0 0 25 0 1 0 224.0 47.2 0.0 0.0 0 0 0 26 0 1 0 139.0 17.0 0.0 0.0 0 0 0 27 0 1 0 281.0 75.5 0.0 0.0 0 0 0 28 0 1 0 206.0 27.6 0.0 0.0 0 0 0 29 0 1 0 283.5 26.9 0.0 0.0 0 0 0 30 2 1.0475 0 0.0 0.0 250.0 0 0 0 0 31 2 0.9820 0 9.2 4.6 0 0 0 0 0 32 2 0.9831 0 0.0 0.0 650.0 0 0 0 0 33 2 0.9972 0 0.0 0.0 632.0 0 0 0 0 34 2 1.0123 0 0.0 0.0 508.0 0 0 0 0 35 2 1.0493 0 0.0 0.0 650.0 0 0 0 0 36 2 1.0635 0 0.0 0.0 560.0 0 0 0 0 37 2 1.0278 0 0.0 0.0 540.0 0 0 0 0 38 2 1.0265 0 0.0 0.0 830.0 0 0 0 0 39 2 1.0300 0 1104.0 250.0 1000.0 0 0 0 0];
linedata=[1 2 0.0035 0.0411 0.6987 0.000 1 39 0.0010 0.0250 0.7500 0.000 2 3 0.0013 0.0151 0.2572 0.000 2 25 0.0070 0.0086 0.1460 0.000 3 4 0.0013 0.0213 0.2214 0.000 3 18 0.0011 0.0133 0.2138 0.000 4 5 0.0008 0.0128 0.1342 0.000 4 14 0.0008 0.0129 0.1382 0.000 5 6 0.0002 0.0026 0.0434 0.000 5 8 0.0008 0.0112 0.1476 0.000 6 7 0.0006 0.0092 0.1130 0.000 6 11 0.0007 0.0082 0.1389 0.000 7 8 0.0004 0.0046 0.0780 0.000 8 9 0.0023 0.0363 0.3804 0.000 9 39 0.0010 0.0250 1.2000 0.000 10 11 0.0004 0.0043 0.0729 0.000 10 13 0.0004 0.0043 0.0729 0.000 13 14 0.0009 0.0101 0.1723 0.000 14 15 0.0018 0.0217 0.3660 0.000 15 16 0.0009 0.0094 0.1710 0.000 16 17 0.0007 0.0089 0.1342 0.000 16 19 0.0016 0.0195 0.3040 0.000 16 21 0.0008 0.0135 0.2548 0.000 16 24 0.0003 0.0059 0.0680 0.000 17 18 0.0007 0.0082 0.1319 0.000 17 27 0.0013 0.0173 0.3216 0.000 21 22 0.0008 0.0140 0.2565 0.000 22 23 0.0006 0.0096 0.1846 0.000 23 24 0.0022 0.0350 0.3610 0.000 25 26 0.0032 0.0323 0.5130 0.000 26 27 0.0014 0.0147 0.2396 0.000 26 28 0.0043 0.0474 0.7802 0.000 26 29 0.0057 0.0625 1.0290 0.000 28 29 0.0014 0.0151 0.2490 0.000 12 11 0.0016 0.0435 0.0000 1.006 12 13 0.0016 0.0435 0.0000 1.006 6 31 0.0000 0.0250 0.0000 1.070 10 32 0.0000 0.0200 0.0000 1.070 19 33 0.0007 0.0142 0.0000 1.070 20 34 0.0009 0.0180 0.0000 1.009 22 35 0.0000 0.0143 0.0000 1.025 23 36 0.0005 0.0272 0.0000 1.000 25 37 0.0006 0.0232 0.0000 1.025 2 30 0.0000 0.0181 0.0000 1.025 29 38 0.0008 0.0156 0.0000 1.025 19 20 0.0007 0.0138 0.0000 1.060];
lfybus % form the bus admittance matrix lfnewton % Power flow solution by Newton-Raphson method busout % Prints the power flow solution on the screen
% Gen. Ra Xd' H gendata=[ 1 0 .006 500 2 0 .0697 30.3 3 0 .0531 35.8 4 0 .0436 28.6 5 0 .132 26 6 0 .05 34.8 7 0 .049 26.4 8 0 .057 24.3 9 0 .057 34.5 10 0 .031 42]; dynamic
but unfortunately, this is a problem in the dimension and I couldn't fix it
this is what I got when I run it :(
Attempted to access xf(20,:); index out of bounds because size(xf)=[15,20].
im really stuck and I need ur help gentelmen

Sign in to answer this question.

Answers (2)

Thomas Koelen
Thomas Koelen on 13 Apr 2015
So your array xf has a size of 15 rows by 20 columns, you are trying to acces row 20, which doesn't excist.
Also, please use the {}code function to make up your post, this is nearly unreadable!
  7 Comments
Show 5 older comments
Sulaiman AlEidan
Sulaiman AlEidan on 20 Apr 2015
I meant if you could try it because im curious to see your modifications
SASWATI MISHRA
SASWATI MISHRA on 17 Nov 2016
In the same code,I get an error like this: Undefined function or method 'ybusbf' for input arguments of type 'double'. Please help me out

Sign in to comment.

rajesh yerupalli
rajesh yerupalli on 1 Jul 2019
Anybody send me simulink model diagram of 10 machine 39 bus system with simulink block names

Categories

Find more on Mathematics and Optimization in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!