电力系统潮流计算 程序及参数NR法，PQ分解法

clear all clc n=14;%请输入节点数; b=20;%请输入支路数; ngnd=1;%请输入对地支路数 tree=xlsread('支路参数'); node1=xlsread('节点参数'); node=node1; gnd=[9 0 0.19];%请输入对地支路（节点号 电导 电纳） nitmax=1000;%最大迭代次数 err=1e-5;%最大允许误差 %对节点类型分类并排序 %1-平衡节点 %2-PV节点 %3-PQ节点 nPQ=0;nPV=0;nBAN=0; %对节点重新排序 exch=zeros(1,n); j=1; for i=1:n if node1(i,2)==3%找出PQ节点个数 nPQ=nPQ+1; node(j,:)=node(i,:);%i,j节点位置互换 exch(i)=j; j=j+1; end end for i=1:n if node1(i,2)==2%找出PV节点个数 nPV=nPV+1; node(j,:)=node1(i,:);%i,j节点位置互换 exch(i)=j; j=j+1; end end for i=1:n if node1(i,2)==1%找出平衡节点 nBAN=i; exch(i)=j; node(n,:)=node1(i,:);%平衡节点放到最后的位置 end end %对应修改支路参数 for j=1:b for m=1:2 tree(j,m)=exch(tree(j,m)); end end %修改对地支路参数 gnd(1)=exch(gnd(1)); %初始化导纳矩阵 YR=zeros(n,n);%电导 YI=zeros(n,n);%电纳 Y=zeros(n,n);%导纳矩阵 %变压器支路的导纳元素 for i=1:b if tree(i,6)==1 mo=tree(i,3)^2+tree(i,4)^2; %自导纳，使用Π型等值（忽略励磁支路） %k*在首端，Z*在末端 YR(tree(i,1),tree(i,1))=YR(tree(i,1),tree(i,1))+(tree(i,3)/mo)/tree(i,7)+(tree(i,3)/mo)*(1-tree(i,7))/tree(i,7)^2; YI(tree(i,1),tree(i,1))=YI(tree(i,1),tree(i,1))+(-tree(i,4)/mo)/tree(i,7)+(-tree(i,4)/mo)*(1-tree(i,7))/tree(i,7)^2; YR(tree(i,2),tree(i,2))=YR(tree(i,2),tree(i,2))+(tree(i,3)/mo)/tree(i,7)+(tree(i,3)/mo)*(tree(i,7)-1)/tree(i,7); YI(tree(i,2),tree(i,2))=YI(tree(i,2),tree(i,2))+(-tree(i,4)/mo)/tree(i,7)+(-tree(i,4)/mo)*(tree(i,7)-1)/tree(i,7); %互导纳 YR(tree(i,1),tree(i,2))=YR(tree(i,1),tree(i,2))-tree(i,3)/(tree(i,7)*mo); YI(tree(i,1),tree(i,2))=YI(tree(i,1),tree(i,2))+tree(i,4)/(tree(i,7)*mo); YR(tree(i,2),tree(i,1))=YR(tree(i,2),tree(i,1))-tree(i,3)/(tree(i,7)*mo); YI(tree(i,2),tree(i,1))=YI(tree(i,2),tree(i,1))+tree(i,4)/(tree(i,7)*mo); end end %线路的导纳元素 for cc=1:b if tree(cc,6)~=1 i=tree(cc,1); j=tree(cc,2); mo=tree(cc,3)^2+tree(cc,4)^2; %自导纳（不计对地电导） YR(i,i)=YR(i,i)+tree(cc,3)/mo; YI(i,i)=YI(i,i)-tree(cc,4)/mo+tree(cc,5); YR(j,j)=YR(j,j)+tree(cc,3)/mo; YI(j,j)=YI(j,j)-tree(cc,4)/mo+tree(cc,5); %互导纳 YR(i,j)=YR(i,j)-tree(cc,3)/mo; YI(i,j)=YI(i,j)+tree(cc,4)/mo; YR(j,i)=YR(j,i)-tree(cc,3)/mo; YI(j,i)=YI(j,i)+tree(cc,4)/mo; end end %添加对地支路 for cc=1:ngnd i=gnd(cc,1); YR(i,i)=YR(i,i)+gnd(cc,2); YI(i,i)=YI(i,i)+gnd(cc,3); end %节点导纳矩阵初始化结束 % 迭代初始化 %计算P-Q给定值 snet=zeros(1,n);%每个节点发电机与负荷的净注入功率 for i=1:n snet(i)=node(i,3)+node(i,4)*1i-node(i,5)-node(i,6)*1i; %功率的给定值 end %（n-1）个P的给定值 %（mPQ）个Q的给定值 f=0;%迭代标志 nit=1;%当前迭代次数 %开始迭代 tic%计时开始 while(f==0&&nit<nitmax) %计算P-Q不平衡量 str=zeros(1,n);%各节点净注入功率 for i=1:n psum=0; qsum=0; for j=1:n %节点电压方程（极坐标形式） psum=psum+node(i,7)*node(j,7)*YR(i,j)*cos(node(i,8)-node(j,8))+node(i,7)*node(j,7)*YI(i,j)*sin(node(i,8)-node(j,8)); %计算每个节点传输有功 qsum=qsum+node(i,7)*node(j,7)*YR(i,j)*sin(node(i,8)-node(j,8))-node(i,7)*node(j,7)*YI(i,j)*cos(node(i,8)-node(j,8)); %计算每个节点传输无功 end str(i)=psum+qsum*1i; end %记录每个节点传输的功率 P=zeros(1,n); Q=zeros(1,n); for i=1:n P(i)=real(str(i));%提取各个节点的有功功率并存入数组 Q(i)=imag(str(i)) ;%提取各个节点的无功功率并存入数组 end DS=snet-str; %得到PQ与给定的偏差 %判断是否满足误差要求 DP=zeros(1,n-1); DQ=zeros(1,nPV); %统计前n-1个节点的有功不平衡量绝对值 for i=1:n-1 DP(i)=abs(real(DS(i))); end %统计PQ节点的无功不平衡量绝对值 for i=1:nPQ DQ(i)=abs(imag(DS(i))); end A=max(DP); B=max(DQ); MAX=max(A,B); if MAX<err f=1;%若误差满足要求，则标志位置1，表示停止迭代 end %判断结束，若不满足，继续迭代 DP=zeros(1,n-1); DQ=zeros(1,nPV); %统计前n-1个节点的有功不平衡量 for i=1:n-1 DP(i)=real(DS(i)); end %统计PQ节点的无功不平衡量 for i=1:nPQ DQ(i)=imag(DS(i)); end %形成PQ不平衡量 delt_PQ=[DP';DQ']; %形成雅可比矩阵 H=zeros(n-1,n-1); L=zeros(nPQ,nPQ); N=zeros(n-1,nPQ); K=zeros(nPQ,n-1); %形成矩阵H for i=1:n-1 for j=1:n-1 if i~=j H(i,j)=-node(i,7)*node(j,7)*(YR(i,j)*sin(node(i,8)-node(j,8))-YI(i,j)*cos(node(i,8)-node(j,8))); else H(i,j)=node(i,7)^2*YI(i,i)+Q(i); end end end %形成矩阵N for i=1:n-1 for j=1:nPQ if i~=j N(i,j)=-node(i,7)*node(j,7)*(YR(i,j)*cos(node(i,8)-node(j,8))+YI(i,j)*sin(node(i,8)-node(j,8))); else N(i,j)=-node(i,7)^2*YR(i,i)-P(i); end end end %形成矩阵K for i=1:nPQ for j=1:n-1 if i~=j K(i,j)=node(i,7)*node(j,7)*(YR(i,j)*cos(node(i,8)-node(j,8))+YI(i,j)*sin(node(i,8)-node(j,8))); else K(i,j)=node(i,7)^2*YR(i,i)-P(i); end end end %形成矩阵L for i=1:nPQ for j=1:nPQ if i~=j L(i,j)=-node(i,7)*node(j,7)*(YR(i,j)*sin(node(i,8)-node(j,8))-YI(i,j)*cos(node(i,8)-node(j,8))); else L(i,j)=node(i,7)^2*YI(i,i)-Q(i); end end end JX=[H N;K L]; %求电压和相角的修正值 Ud2=zeros(nPQ,nPQ); for i=1:nPQ for j=1:nPQ if i==j Ud2(i,j)=node(i,7); end end end delt=-inv(JX)*delt_PQ; %分别求幅值和相位的修正量 delt_p=delt(1:n-1,1); delt_a=Ud2*delt(n:n+nPQ-1,1); %求得电压和相位的修正量 node(1:n-1,8)=node(1:n-1,8)+delt_p; node(1:nPQ,7)=node(1:nPQ,7)+delt_a; nit=nit+1; end %迭代结束，得到每个节点电压幅值和相位的结果 toc%计时结束 %开始计算发电机功率 %负荷功率均为给定值 %str中记录了结果中每个节点注入的功率 for i=1:n if node(i,2)==2%PV节点，需求解注入的无功 node(i,4)=imag(str(i))+node(i,6); elseif node(i,2)==1%平衡节点，需求解注入的有功无功 node(i,3)=real(str(i))+node(i,5); node(i,4)=imag(str(i))+node(i,6); end end gnd tree node=sortrows(node) xlswrite('计算结果.xlsx',node,'A2:H15');