基于Matlab实现含弱环网的配电网潮流计算

% function [objectfun,loss_sum] = runpf_test33(Z,par) %该程序用前推回推法计算配电网潮流% clc clear par=[5 6 8];%设定并联电容器组的初始投入组数 Urat=1.0;%额定电压幅值 % myresult=fopen('D:\matlab2009\work\powerflow\output.txt','wt'); Xm=3.54708;X1=0.20891;%假设的风机参数值 P=300/600;%在给定风速下的风机功率(标么值) Uf=0.42/0.69;%风机给定电压 % c=1.07;%定义风机迭代误差 delat=sqrt(Uf^4-4*(P^2)*(X1^2)); if delat>=0 Q=Uf^2/Xm+(Uf^2-delat)/(2*X1);%计算风机无功功率 else disp('错误'); end a1=P/sqrt(P^2+Q^2);%未投入并联电容器组时的功率因数 %a1=0.89;%风机的额定功率因数，未加电容器组时 %P_angle1=acos(a1); Q_Unit=40;%每台并联电容器的容量（有名值） U_N=1.0;%电容器额定工作电压 V1=0.639/0.69;%电容器初始电压 a2=0.92;%要求达到的功率因数 Qnn1=(par(1))*(Q_Unit)*(V1/U_N)^2; Qnn2=(par(2))*(Q_Unit)*(V1/U_N)^2; Qnn3=(par(3))*(Q_Unit)*(V1/U_N)^2; % P_angle2=acos(a2); % Qc=P/tan(P_angle2); % %Qc=Pe*(sqrt(1/(a1)^2-1)-sqrt(1/(a2)^2-1));%并联电容器输出的无功补偿量 % Q_Unit=40;%每台并联电容器的容量（有名值） % N_needf=round(((Qc-Q)*600)/Q_Unit);%需要的并联电容器组数,取整 % Qn=(N_needf)*(Q_Unit)*(V1/U_N)^2;%Ud应该如何循环调用？？？？？？？？？？？？？ % a3=P/sqrt(P^2+(Qc-Q)^2);%带并联电容器组的风机的功率因数 P=-P*600; %Q=Q*10000; Q1f1=-Q*600+Qnn1; Q1f2=-Q*600+Qnn2; Q1f3=-Q*600+Qnn3; %Num=input('输入风机插入处节点号Num='); Num =21; Num1=25; Num2=17; Num3=21;%并联电容器节点号 Num4=25; Num5=17; n =33; b = 32; Sb = 10; Ub = 12.66; %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; G=Z(Num,5);G1=Z(Num1,5);G2=Z(Num2,5);%风机接入节点数 Z(Num,5)=(real(G)+P)+i*(imag(G)-Q1f1); Z(Num1,5)=(real(G1)+P)+i*(imag(G1)-Q1f2); Z(Num2,5)=(real(G2)+P)+i*(imag(G2)-Q1f3); M=Z(Num3,5);M1=Z(Num4,5);M2=Z(Num5,5); %Z(Num,5)=P+60+i*(Q1+25); Z(:,5)=Z(:,5)/Sb/1000;%handle datas V=ones(n,1); %%%%循环变量初始化%%%% Q1=zeros(1,3); delat1=zeros(1,3); N_need=zeros(1,3); Qn1=zeros(1,3); Qf=zeros(1,3); U1=zeros(n,1); U2=zeros(n,1); R=zeros(1,3);%带并联电容器组的风电机组的初始功率因数 c=zeros(1,3);%风电机组的初始功率因数 %迭代开始处 t=0; k=0; while t<b & k<20 %计算节点注入电流 for L=1:b a=Z(L,3); ua=V(a+1,1); I(a,1)=conj(Z(a,5)/ua); end % I(Num3,1)=I(Num3,1)+U2(Num3,1)/R; % I(Num4,1)=I(Num4,1)+U2(Num4,1)/R; % I(Num5,1)=I(Num5,1)+U2(Num5,1)/R; %回推算支路电流 J=zeros(b,1); % Ploss=zeros(b,1);%有功损耗初始化 % Qloss=zeros(b,1);%无功损耗初始化 % U=zeros(k+1,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 a=Z(L,3)+1; a1=Z(L,2)+1; V(a,1)=V(a1,1)-Z(L,4)*J(L,1); % U(a,1)=U(a-1,1)+1*( V(a,1)-U(a-1,1));%加速迭代的速度 % V(a,1)=U(a,1); end U1=V(:,1); %收敛判定 t=0; for a=2:n SS=V(a,1)*conj(I(a-1,1)); dp=real(SS-Z(a-1,5)); dq=imag(SS-Z(a-1,5)); S(a-1,1)=SS; ddp=abs(dp); ddq=abs(dq); L1=(ddp<0.0000001)&(ddq<0.0000001); F1(a-1,1)=L1; if L1==1 t=t+1; end end k=k+1; Unum=[abs(V(Num,1)) abs(V(Num1,1)) abs(V(Num2,1))]; P=300/600;%在给定风速下的风机功率(标么值) for r=1:1:3 delat1(r)=sqrt(Unum(r)^4-4*(P^2)*(X1^2)); if delat1(r)>=0 Qf(r)=Unum(r)^2/Xm+(Unum(r)^2-delat1(r))/(2*X1);%计算风机无功功率 else disp('错误'); end % N_need(r)=abs(round(((Qc-Qf(r))*600)/Q_Unit));%需要的并联电容器组数,取整 % Ur=0.639;%电容器初始电压 aa1=(P*600)/sqrt((P*600)^2+(par(1)*40-Qf(1)*600)^2);%带并联电容器组的风机的功率因数 aa2=P/sqrt(P^2+(par(2)*40-Qf(2))^2); aa3=P/sqrt(P^2+(par(3)*40-Qf(3))^2); R=[aa1 aa2 aa3]; if R(1)>0.9&R(2)>0.9&R(3)>0.9 par=[par(1)-1 par(2)-1 par(3)-1]; else par=[par(1)+1 par(2)+1 par(3)+1]; end Qn1(r)=(par(r))*(Q_Unit)*(Unum(r)/U_N)^2;%Ud应该如何循环调用？？？？？？？？？？？？？ % a3=P/sqrt(P^2+(Qc-Q)^2);%带并联电容器组的风机的功率因数 Q1(r)=-Qf(r)*600+Qn1(r); c(r)=(P*600)/sqrt((P*600)^2+Q1(r)^2); end P=-P*600; % % Z(Num,5)=(real(Z(Num,5))+P)+i*(imag(Z(Num,5))+Q1); %%%%%并联风机%%%%%%% Z(Num,5)=(real(G)+P)+i*(imag(G)-Q1(1)); Z(Num1,5)=(real(G1)+P)+i*(imag(G1)-Q1(2)); Z(Num2,5)=(real(G2)+P)+i*(imag(G2)-Q1(3)); Z(Num,5)=Z(Num,5)/Sb/1000;%化为标幺值 Z(Num1,5)=Z(Num1,5)/Sb/1000;%化为标幺值 Z(Num2,5)=Z(Num2,5)/Sb/1000;%化为标幺值 % par=discrete33N(par); %%%并联电容器%%%% % disp(par) % Z(Num3,5)=real(M)+i*(imag(M)-par(1)*10); % Z(Num4,5)=real(M1)+i*(imag(M1)-par(2)*10); % Z(Num5,5)=real(M2)+i*(imag(M2)-par(3)*10); % Z(Num3,5)=Z(Num3,5)/Sb/1000;%化为标幺值 % Z(Num4,5)=Z(Num4,5)/Sb/1000;%化为标幺值 % Z(Num5,5)=Z(Num5,5)/Sb/1000;%化为标幺值 U2=V(:,1);%保留第k-1次迭代数值 end % t2=clock; % T=t2-t1; % costime=T; %输出结果与否 % disp('输出直角坐标各节点电压'); % disp(V); disp('显示迭代次数'); disp(k); % disp('显示收敛节点情况,"1"表示收敛，"0"表示不收敛'); % disp(F1); %disp(U(a,1)); for a=1:b if F1(a,1)==0 disp('显示不收敛节点号、计算功率'); disp(a);disp(S(a,1)); end end for a=1:n Vm(a,1)=abs(V(a,1));Va(a,1)=angle(V(a,1)); end disp('输出各节点电压幅值'); disp(Vm); % disp('输出各节点电压相角'); % disp(Va) for Line=1:b Ploss(Line,1)=abs((Z(Line,4))*((real(Z(Line,5)))^2+(imag(Z(Line,5)))^2))/(Vm(Line,1)^2); % Ploss(Line,1)=3*real(Z(Line,4)*(J(Line,1).*conj(J(Line,1))));%有功损耗 end loss_sum=sum(Ploss(:,1)); % Loss2=sum(Ploss(:,1)); % disp('loss_sum数值') % disp(loss_sum); U=abs(Vm); objectfun=loss_sum+500*max(abs(U-Urat)); % disp('objectfun数值1') % disp(objectfun) % end