IEEE33节点的并网燃料电池发电.zip

%该程序用前推回推法计算配电网潮流% clear clc % global myresult Q P Num k a b % myresult=fopen('D:\matlab2009\work\powerflow\output.txt','wt'); % Xc=4.105;Xm=3.54708;V=10;%假设的风机参数值 % X=0.20891;P=1500; % b=1.07;%定义风机迭代误差 % Q=V^2*((Xc-Xm)/(Xc*Xm))+(X/V^2)*P^2;%计算风机无功功率 % Num=input('输入风机插入处节点号Num='); % n=input('请输入节点数n='); %%%%%%%%并网燃料电池参数及其初始化%%%%%%%% Pn=150/10000;%有功功率输出恒定值（p.u） Qmin=0;%燃料电池无功功率的下限 Rmin=0.75;%最小功率因数 Smax=Pn/Rmin;%变换器的最大容量 Qmax=sqrt(Smax^2-Pn^2);%燃料电池的无功输出上限 % Q1=(Qmax+Qmin)/2;%燃料电池并网节点无功功率初始化 Q1=0;%燃料电池初始 Q2=0; Q3=0; %%%%%%%%%%%%设定燃料电池节点电压值%%%%%%% Vs=0.995;%（标幺值） %%%%%%%%%%燃料电池并网节点%%%%%%%%%%%%% Num2=23; Num3=29; %%%%%%%%%系统参数%%%%%%%%%%%%% n=33; % b=input('请输入支路数b='); b=32; %disp('请输入支路阻抗和节点功率矩阵，第一列存支路号，第二列存首节点号，第三列存尾节点号，第四列存支路自阻抗,第五列存尾节点给定功率'); % Sb=input('请输入三相功率基准值sb='); Sb=10; % Ub=input('请输入线电压的基准值Ub='); Ub=12.66;%（kV) % Z=input('请输入母线系统支路参数和负荷数据形成的矩阵Z='); Z = [1 0 1 0.0922+i*0.0470 100+i*60 2 1 2 0.4930+i*0.2511 90+i*40 3 2 3 0.3660+i*0.1864 120+i*80 4 3 4 0.3811+i*0.1941 60+i*30 5 4 5 0.8190+i*0.7070 60+i*20 6 5 6 0.1872+i*0.6188 200+i*100 7 6 7 0.7114+i*0.2351 200+i*100 8 7 8 1.0300+i*0.7400 60+i*20 9 8 9 1.0440+i*0.7400 60+i*20 10 9 10 0.1966+i*0.0650 45+i*30 11 10 11 0.3744+i*0.1238 60+i*35 12 11 12 1.468+i*1.1550 60+i*35 13 12 13 0.5416+i*0.7129 120+i*80 14 13 14 0.5910+i*0.5260 60+i*10 15 14 15 0.7463+i*0.5450 60+i*20 16 15 16 1.2890+i*1.7210 60+i*20 17 16 17 0.7320+i*0.5740 90+i*40 18 1 18 0.1640+i*0.1565 90+i*40 19 18 19 1.5042+i*1.3554 90+i*40 20 19 20 0.4095+i*0.4784 90+i*40 21 20 21 0.7089+i*0.9373 90+i*40 22 2 22 0.4512+i*0.3083 90+i*50 23 22 23 0.8980+i*0.7091 420+i*200 24 23 24 0.8960+i*0.7011 420+i*200 25 5 25 0.2030+i*0.1034 60+i*25 26 25 26 0.2842+i*0.1447 60+i*25 27 26 27 1.0590+i*0.9337 60+i*20 28 27 28 0.8042+i*0.7006 120+i*70 29 28 29 0.5075+i*0.2585 200+i*600 30 29 30 0.9744+i*0.9630 150+i*70 31 30 31 0.3105+i*0.3619 210+i*100 32 31 32 0.3410+i*0.5302 60+i*40]; Zb=Ub^2/Sb; Z(:,4)=Z(:,4)/Zb; Y=Z(Num2,5);%燃料电池接入节点数 Y1=Z(Num3,5); Z(Num2,5)=(real(Y)-Pn*10000)+i*(imag(Y)-Q1*10000); Z(Num3,5)=(real(Y1)-Pn*10000)+i*(imag(Y1)-Q1*10000); Z(:,5)=Z(:,5)/Sb/1000;%handle datas % Z(Num,5)=(P+i*Q)/1000; k=0; V=ones(n,1); t=0; Iinject=0; Iinject1=0; %迭代开始处 while t<b & k<100 %计算节点注入电流 for l=1:b j=Z(l,3); ua=V(j+1,1); I(j,1)=conj(Z(j,5)/ua); end I(Num2,1)=conj(Z(Num2,5)/V(Num2+1,1)+i*Iinject); I(Num3,1)=conj(Z(Num3,5)/V(Num3+1,1)+i*Iinject1); % for l=Num2+1:b % j=Z(l,3); % ua=V(j+1,1); % I(j,1)=conj(Z(j,5)/ua); % end %回推算支路电流 J=zeros(b,1); U=zeros(k,1); l=b; J(l)=J(l)+I(l); for jj=1:b-1 l=l-1; for m=l+1:b if Z(m,2)==Z(l,3) J(l)=J(l)+J(m); end end J(l)=J(l)+I(l); end %前推算节点电压 for l=1:b j=Z(l,3)+1; a=Z(l,2)+1; V(j,1)=V(a,1)-Z(l,4)*J(l,1); end Vf=abs(V(:,1)); %收敛判定 t=0; for j=2:n SS=V(j,1)*conj(I(j-1,1)); dp=real(SS-Z(j-1,5)); dq=imag(SS-Z(j-1,5)); S(j-1,1)=SS; ddp=abs(dp); ddq=abs(dq); L1=(ddp<0.0000001)&(ddq<0.0000001); F(j-1,1)=L1; if L1==1 t=t+1; end end k=k+1; if Q2>=Qmax Iinject=0; Q2=Qmax; Z(Num2,5)=(real(Y)-Pn*10000)+i*(imag(Y)-Q2*10000); Z(Num2,5)=Z(Num2,5)/Sb/1000;%化为标幺值 else Vd=abs(Vf(Num2,1)-Vs); Iinject=Vd/imag(Z(Num2,4)); Qd=imag(V(Num2,1)*conj(I(Num2,1))); Q2=imag(Y/10000)+Qd; end if Q3>=Qmax Iinject1=0; Q3=Qmax; Z(Num3,5)=(real(Y1)-Pn*10000)+i*(imag(Y1)-Q3*10000); Z(Num3,5)=Z(Num3,5)/Sb/1000;%化为标幺值 else Vd1=abs(Vf(Num3,1)-Vs); Iinject1=Vd1/imag(Z(Num3,4)); Qd1=imag(V(Num3,1)*conj(I(Num3,1))); Q3=imag(Y1/10000)+Qd1; end end L=zeros(32,1); L(1,1)=4*((imag(Z(1,4))*real(Z(1,5))-real(Z(1,4))*imag(Z(1,5)))^2+(imag(Z(1,4))*imag(Z(1,5))+real(Z(1,4))*real(Z(1,5)))); for nn=2:b jj=Z(nn,2)+1; L(nn,1)=(4*((imag(Z(nn,4))*real(Z(nn,5))-real(Z(nn,4))*imag(Z(nn,5)))^2+(imag(Z(nn,4))*imag(Z(nn,5))+real(Z(nn,4))*real(Z(nn,5)))*((abs(V(jj,1))))^2))/(abs(V(jj,1)))^4; Lmax=max(L(:,1)); end %输出结果与否 % disp('输出直角坐标各节点电压'); % disp(V); disp('显示迭代次数'); disp(k); disp('显示收敛节点情况,"1"表示收敛，"0"表示不收敛'); disp(F); % disp(U(a)); for j=1:b if F(j,1)==0 disp('显示不收敛节点号、计算功率'); disp(j);disp(S(j,1)); end end for j=1:n Vm(j,1)=abs(V(j,1));Va(j,1)=angle(V(j,1)); end disp('输出各节点电压幅值'); disp(Vm); disp('输出各节点电压相角'); disp(Va) for Line=1:b % Ploss(Line,1)=real(Z(Line,4))*(Z(Line,5).*conj(Z(Line,5))/(Vm(Line,1)^2)); Ploss(Line,1)=3*real(Z(Line,4)*(J(Line,1).*conj(J(Line,1))));%有功损耗 end sunhao1=sum(Ploss(:,1));%总和 %返回开始处