I'm trying to solve this but it says ( Attempted to access x(2); index out of bounds because numel(x)=1 ) ..I just want to change the solver from ode32 '' I marked it as a comment'' to modified euler method but it showing this problem
1 view (last 30 days)
Show older comments
global Pm f H E Y th ngg f=60; R=linedata(:,3); X=linedata(:,4); % added 18/10/2012 %zdd=gendata(:,2)+j*gendata(:,3); ngr=gendata(:,1); %H=gendata(:,4); ngg=length(gendata(:,1)); %% for k=1:ngg zdd(ngr(k))=gendata(k, 2)+j*gendata(k,3); %H(ngr(k))=gendata(k, 4); H(k)=gendata(k,4); % new end %% for k=1:ngg I=conj(S(ngr(k)))/conj(V(ngr(k))); %Ep(ngr(k)) = V(ngr(k))+zdd(ngr(k))*I; %Pm(ngr(k))=real(S(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) = gendata(k, 2); %X(nbr+k) = gendata(k, 3);
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. -> '); 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 %%%%%%% from here the changing of solver from ode32 to meuler_simple.%%%%%%% 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('dfpek', tspan, x0);
d= 0.02;N=20; t1=0; h=(tc-t1)/N; for k=1:N t1=t1+h; [t1,x]= meuler_simple(@dfpek, t1,x0,tc,N); x1=x; x1=[x1; x0]; end
x1c =x1(length(x1), :);
....................................... Y=abs(Yaf); th=angle(Yaf); ....................................... % tspan = [tc, tf]; % [t2,xc] =ode23('afpek', tspan, x0c);
l=(tf-tc)/N; t2=tc for n=1:N t2=t2+l; [t2,x]= meuler_simple(@afpek, tc,x1c,tf,N); x2=x; x2=[x2; x1c] end
t =[t1; t2]; x = [x1; x2];
0 Comments
Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!