PQ分解法计算大电网潮流程序

functionPQ%用PQ分解法计算大电网潮流%%bus数组1.节点编号2.节点电压3.节点电压角度4.注入有功5.注入无功6.节点类型(1PQ2PV3平衡)%line数组1.始端节点编号2.末端节点编号3.电阻4电抗5电导G6电纳B7.变比%打开数据文件clearclcbus=load('ieee14bus.txt');line=load('ieee14line.txt');linenum(:,[1,2])=line(:,[1,2]);[nb,~]=size(bus);[nl,~]=size(line);%nodenum=[(1:nb)'bus(:,1)];%带入子函数数据处理[bus,line,nPQ,nPV,nSW,nodenum]=change1_busline(bus,line);%对节点重新编号Y=admittance(bus,line,1);%生成节点导纳矩阵Y1=admittance(bus,line,2);%生成化简条件3的矩阵B1Y2=admittance(bus,line,3);%生成化简条件3的矩阵B2%-----------------------------------------------------%%临时添加的测试数据%nPQ=4;nPV=0;nSW=1;nb=5;%Y=[10.834-32.5i-1.667+5i-1.667+5i-2.5+7.5i-5+15i%-1.667+5i12.917-38.75i-10+30i0-1.25+3.75i%-1.667+5i-10+30i12.917-38.75i-1.25+3.75i0%-2.5+7.5i0-1.25+3.75i3.75-11.25i0%-5+15i-1.25+3.75i006.25-18.75i];%%Y1=[10.834-32.5i-1.667+5i-1.667+5i-2.5+7.5i-5+15i%-1.667+5i12.917-38.75i-10+30i0-1.25+3.75i%-1.667+5i-10+30i12.917-38.75i-1.25+3.75i0%-2.5+7.5i0-1.25+3.75i3.75-11.25i0%-5+15i-1.25+3.75i006.25-18.75i];%%bus=[1100.20.21%210-0.45-0.151%310-0.4-0.051%410-0.6-0.11%51.060003];%%line=[521.25-3.75000%2310-30000%341.25-3.75000%412.5-7.5000%121.667-5000%131.667-5000%155-15000];%-------------------------------------------------------bus_PV0=bus((nPQ+1):end,2)';%1.05*ones(1,nPV+nSW);bus_U=[ones(1,nPQ)bus_PV0]';%电压幅值bus_e=zeros(nb,1);%电压角度delta_P=zeros(nPQ+nPV,1);delta_Q=zeros(nPQ,1);%delta_e=zeros(nb-1,1);%delta_U=zeros(nPQ,1);c=0;KP=1;KQ=1;%KPKQ用来判断有功、无功是否收敛G=real(Y);B=imag(Y);B10=imag(Y1);B20=imag(Y2);%矩阵B0是进行化简三后的节点导纳矩阵虚部%形成解耦潮流的系数矩阵B1和B2B1=B10(1:nb-1,1:nb-1);B2=B20(1:nPQ,1:nPQ);whilec80%求解PQ的不平衡量forii=1:nPQ+nPVdelta_P(ii)=bus(ii,4);forjj=1:nbdelta_P(ii)=delta_P(ii)-bus_U(ii)*bus_U(jj)*(G(ii,jj)*cos(bus_e(ii)-bus_e(jj))+B(ii,jj)*sin(bus_e(ii)-bus_e(jj)));endendUP=diag(bus_U(1:(nb-1)));%U矩阵利用各节点电压形成对角阵，来计算修正方程，对角线上的元素与bus_U列元素一一对应error_P=UP\delta_P;ifmax(abs(error_P))0.00001delta_e=-(UP*B1)\error_P;bus_e=bus_e+[delta_e;0];c=c+1;KQ=1;elseKP=0;ifKQ~=0elsebreakendendforii=1:nPQdelta_Q(ii)=bus(ii,5);forjj=1:nbdelta_Q(ii)=delta_Q(ii)-bus_U(ii)*bus_U(jj)*(G(ii,jj)*sin(bus_e(ii)-bus_e(jj))-B(ii,jj)*cos(bus_e(ii)-bus_e(jj)));endendUQ=diag(bus_U(1:(nb-nPV-nSW)));error_Q=UQ\delta_Q;ifmax(abs(error_Q))0.00001delta_U=-B2\error_Q;bus_U=bus_U+[delta_U;zeros((nPV+nSW),1)];c=c+1;KP=1;elseKQ=0;ifKP~=0elsebreakendendend%至此得到收敛的节点电压值%----------------------------------------------%--------------------------------------------%对计算结果进行数据处理%将节点结果用原节点编号表示bus_Ue=zeros(nb,3);bus_Ue(:,[1,2,3])=[nodenum(:,2)bus_Ubus_e/pi*180];forii=1:nbforjj=ii+1:nbifbus_Ue(ii,1)bus_Ue(jj,1)t=bus_Ue(ii,:);bus_Ue(ii,:)=bus_Ue(jj,:);bus_Ue(jj,:)=t;endendend%r_U是收敛的电压表达成复数的形式r_U=zeros(nb,1);fork=1:nbr_U(k)=bus_U(k)*(cos(bus_e(k))+1i*sin(bus_e(k)));end%计算平衡节点功率SW_S=0;SW_S=SW_S+r_U(nb)*conj(Y(nb,:))*conj(r_U);%计算各支路功率Sijline_S=zeros(nb,nb);line_S0=zeros(nb,nb);forii=1:nbforjj=1:nbline_S(ii,jj)=r_U(ii)*(conj(r_U(ii))*conj(Y(ii,ii))+(conj(r_U(ii))-conj(r_U(jj)))*conj(Y(ii,jj)));endend%------------------------------------------%把线路结果还原成原节点编号对应的结果forii=1:nbforjj=1:nbline_S0(nodenum(ii,2),nodenum(jj,2))=line_S(ii,jj);endendline_P=real(line_S0);line_Q=imag(line_S0);%计算各支路损耗delta_S=zeros(nl,1);fork=1:nla=linenum(k,1);b=linenum(k,2);delta_S(k)=line_S0(a,b)+line_S(b,a);end%计算网络总损耗S0=sum(delta_S);%-------------------------------------------------%将计算结果输入指定文件fid=fopen('C:\Users\lr\Desktop\matlab练习\训练题\大电网潮流计算\ieee14_out.txt','wt');fprintf(fid,'节点号\t节点电压幅值\t节点电压角度\n');fork=1:nbfprintf(fid,'%d\t%f\t%f\n',k,bus_Ue(k,1),bus_Ue(k,2));endfprintf(fid,'支路首端\t支路末端\t支路有功\t支路无功\t支路损耗\n');fork=1:nlfprintf(fid,'%d\t\t%d\t\t%f\t%f\t%f\n',linenum(k,1),linenum(k,2),line_P(linenum(k,1),linenum(k,2)),line_Q(linenum(k,1),linenum(k,2)),delta_S(k));endfprintf(fid,'平衡节点功率=%f\n',SW_S);fprintf(fid,'网络总损耗=%f\n',S0);fclose(fid);endfunction[bus,line,nPQ,nPV,nSW,nodenum]=change1_busline(bus,line)%此函数用来对原始输入节点、线路数据进行重新编号%%bus数组1.节点编号2.节点电压3.节点电压角度4.注入有功5.注入无功6.节点类型(1PQ2PV3平衡)%line数组1.始端节点编号2.末端节点编号3.电阻4电抗5电导G6电纳B7.变比[nb,~]=size(bus);[nl,~]=size(line);%nodenum=[(1:nb)'bus(:,1)];nPQ=0;nPV=0;nSW=0;%PQ=[];PQ=zeros(nb,6);PV=zeros(nb,6);SW=zeros(nb,6);%PQPV平衡节点的个数fork=1:nbswitchbus(k,6)case1nPQ=nPQ+1;PQ(nPQ,:)=bus(k,:);case2nPV=nPV+1;PV(nPV,:)=bus(k,:);case3nSW=nSW+1;SW(nSW,:)=bus(k,:);otherwisedisp('节点数据类型出错！');endend%生成重新编号后的节点数据矩阵bus=[PQ;PV;SW];nodenum=[(1:nb)'bus(:,1)];%第一列为新的节点编号，第二列为对应的旧节点编号bus(:,1)=(1:nb)';%至此实现了节点数据的重新编号%------------------------------------------------------%对线路数据重新编号%nodenum=[(1:nb)'bus(:,1)];%第一列为新的节点编号，第二列为对应的旧节点编号forii=1:nl[r1,~]=find(nodenum(:,2)==line(ii,1));line(ii,1)=nodenum(r1,1);[r2,~]=find(nodenum(:,2)==line(ii,2));line(ii,2)=nodenum(r2,1);endendfunctionY=admittance(bus,line,c)%此函数用来形成节点导纳矩阵%%bus数组1.节点编号2.节点电压3.节点电压角度4.注入有功5.注入无功6.节点类型(1PQ2PV3平衡)%line数组1.始端节点编号2.末端节点编号3.电阻4电抗5电导6电纳B/27.变比%c是用来控制形成节点导纳矩阵的方式的，c=1形成一般的节点导纳矩阵，可以用来确定B2，c=2,形成化简条件3的节点导纳矩阵，确定B1[nb,~]=size(bus);[nl,~]=size(line);Y=zeros(nb,nb);zt=zeros(nl,1);yt=zeros(nl,1);ym=zeros(nl,1);I=zeros(nl,1);J=zeros(nl,1);K=zeros(nl,1);switchccase1fork=1:nlzt(k)=line(k,3)+1