coding new england 39-bus system

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

Answer 1

[Thomas Koelen]

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!

[Sulaiman AlEidan]

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

| | | | | | this is the data for the 39-bus system | | | | | | *

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

first of all Im sorry because I wanted to arrange it but i didnt know about {code}

so, what do you suggest to change ?

thanks a lot for your reply

[Thomas Koelen]

At what line does the error occur?

[Sulaiman AlEidan]

line 92 which is this

x0c =xf(length(xf), :); Y=abs(Yaf); th=angle(Yaf); %[t2,xc] =ode23('afpek', tc, tf, x0c, tol); % Postfault solution (use with MATLAB 4)

[Thomas Koelen]

create a variable that's size(xf) then change length(xf) to variable(1,1) or (1,2), try both and check if it works.

[Sulaiman AlEidan]

I tried them both and it didn't work ,, :(

can you please run it ?

Maybe you meant something else.

I really appreciate your help ,,

thanks a lot

[Sulaiman AlEidan]

I meant if you could try it because im curious to see your modifications

[SASWATI MISHRA]

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

Answer 2

[rajesh yerupalli]

Anybody send me simulink model diagram of 10 machine 39 bus system with simulink block names