Decouple.m
资源名称:power.rar [点击查看]
上传用户:eighthdate
上传日期:2014-05-24
资源大小:270k
文件大小:4k
源码类别:
其他行业
开发平台:
Matlab
- % Fast Decoupled Power Flow Solution
- % Copyright (c)1998 by Hadi Saadat.
- ns=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
- 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;
- % Start of iterations
- while maxerror >= accuracy & iter <= maxiter % Test for max. power mismatch
- iter = iter+1;
- id=0; iv=0;
- for n=1:nbus
- nn=n-nss(n);
- J11=0; J33=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));
- 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) - Qsh(n);
- if Qmax(n) ~= 0
- if iter <= 20 % Between the 1th & 6th iterations
- if iter >= 10 % the Mvar of generator buses are
- if Qgc < Qmin(n), % tested. If not within limits Vm(n)
- Vm(n) = Vm(n) + 0.005; % is changed in steps of 0.05 pu to
- elseif Qgc > Qmax(n), % bring the generator Mvar within
- Vm(n) = Vm(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);
- elseif kb(n) ==2
- S(n)=P(n)+j*Q(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 Pk Qk DP DQ J11 J33 B1 B1inv B2 B2inv DPV DQV Dd delta ib id ii iv jb jj