目录

引言求解过程（1）求解主函数（2）雅克比矩阵求解函数（3）线损计算程序（4）三相短路的短路容量计算程序

引言

“N－1”准则 电网供电安全所采用的准则，要求同级电压网络正常运行方式下任一元件（变压器、线路、母线）无故障或因故障断开时，电网应能稳定运行和不损失负荷正常供电，其它元件不过负荷，电压和频率均在允许范围内。

求解过程

（1）求解主函数

% clc % V_fangshi1=[1.05,1.0207,1.0384,1.0254,1.0306,1.0455]; % for Ls=1:8 % 线路参数的有名值计算------------------------------------------------------ % 正序参数 L=[20,21,17,25,27,30,30,30]; R=0.054*L; X=0.308*L; B_2=pi*50*0.0116*10^-6*L; R_cocy=R; B_copy=B_2; R(4:6)=R(4:6)/2; X(4:6)=X(4:6)/2; B_2(4:6)=B_2(4:6)*2; Canshu_zhengxu=[R;X;B_2]; % 零序参数 R_0=0.204*L; X_0=0.968*L; B_2_0=pi*50*0.0078*10^-6*L; R_0(4:6)=R_0(4:6)/2; X_0(4:6)=X_0(4:6)/2; B_2_0(4:6)=B_2_0(4:6)*2; Canshu_0_zhengxu=[R_0;X_0;B_2_0]; % 线路参数的标幺值计算 R_B=R*(100/(220*220)); X_B=X*(100/(220*220)); B_2_B=B_2*((220*220)/100); R_0_B=R_0*(100/(220*220)); X_0_B=X_0*(100/(220*220)); B_2_0_B=B_2_0*((220*220)/100); Canshu_biaoyazhi_lingxu=[R_B;X_B;B_2_B]; Canshu_0_biaoyazhi_lingxu=[R_0_B;X_0_B;B_2_0_B]; Canshu=[Canshu_zhengxu;Canshu_0_zhengxu;Canshu_biaoyazhi_lingxu;Canshu_0_biaoyazhi_lingxu]; % 潮流计算主程序------------------------------------------------------------- %程序解释 %节点类型 标号 %PQ节点 2,3,4,5 %PV节点 6 %平衡节点 1 %能计算给定标幺值网络，有且仅有一个平衡节点的潮流 %注意：母线标号顺序要求：PQ节点-PV节点-平衡节点 %若某元件不存在，其导纳为0，阻抗为inf % 第一步 按已知参数形成节点导纳矩阵------------------------------------------- % nl=input('请输入节点支路数'); %nl是支路数； % n=input('请输入节点数'); %n是节点数； nl=8; %nl是支路数； n=6; %n是节点数； %--------------------------------------------------------------------------- % Ls=input('请输入N-1检验的支路数'); %注：这是N-1程序，平时注释掉，需要时取消注释； % R_B(1)=2*R_B(1); %1:1-4,2 % X_B(1)=2*X_B(1); % B_2_B(1)=B_2_B(1)/2; % R_B(3)=2*R_B(3); %1:1-4,2 % X_B(3)=2*X_B(3); % B_2_B(3)=B_2_B(3)/2; %--------------------------------------------------------------------------- b=complex(0,B_2_B); %b是支路导纳； Z=complex(R_B,X_B); %Z是支路阻抗； Y=zeros(6,6); %Y是节点导纳矩阵 B1=[1,4,1,1,6,2,2,5; %这个矩阵是有支路两端的号数得到的，这个矩阵的建立能简化节点导纳矩阵的形成； 4,5,5,3,3,6,4,6]; B1_copy=B1; %-------------------------------------------------------------------------- % Ls=input('请输入N-1检验的支路数'); % B1(:,Ls)=[]; % Z(:,Ls)=[]; % b(:,Ls)=[]; %------------------------------------------------------------------------ B1=[B1;Z;b]; for i=1:nl p=B1(1,i); q=B1(2,i); Y(p,q)=-1./B1(3,i); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1./B1(3,i)+B1(4,i); Y(q,q)=Y(q,q)+1./B1(3,i)+B1(4,i); end % 第二步 给定节点电压初值 e=zeros(1,n); f=zeros(1,n); for i=1:n e(i)=1.0; f(i)=0; end e(1)=1.05; %PV节点 e(6)=23/22; %平衡节点 % 第三步 计算delta Pi，delta Ui，delta Ui2 G=real(Y); B=imag(Y); %-------------------------------------------------------------------------- Ps=[7.65,-4,-3.5,-4.2,-4.2]; %方式1 Qs=[0,-2.2,-1.8,-1.4,-1.5]; % Ps=[2.55,-2,-1.75,-2.1,-2.1]; % Qs=[0,-1.1,-0.9,-0.7,-0.75]; %-------------------------------------------------------------------------- S1=0; S2=0; DP=zeros(1,n-1); DQ=zeros(1,n-1); % De=zeros(1,n); % Df=zeros(1,n); DW=zeros(2*(n-1),1); DV=zeros(2*(n-1),1); E_record=[]; F_record=[]; % 进入迭代过程-------------------------------------------------------------- ac=0.0000001; %ac是误差； t=0; %t是迭代次数； % PQ节点的计算 while 1 for i=2:n-1 for j=1:n S1=S1+G(i,j)*e(j)-B(i,j)*f(j); S2=S2+G(i,j)*f(j)+B(i,j)*e(j); end DP(i)=Ps(i)-(e(i)*S1+f(i)*S2); DQ(i)=Qs(i)-(f(i)*S1-e(i)*S2); S1=0; S2=0; end % PV节点的计算 for j=1:n S1=S1+G(1,j)*e(j)-B(1,j)*f(j); S2=S2+G(1,j)*f(j)+B(1,j)*e(j); end DP1=Ps(1)-(e(1)*S1+f(1)*S2); DV2=1.05*1.05-e(1)*e(1)-f(1)*f(1); S1=0; S2=0; % 第四步 计算雅克比矩阵 for i=2:n-1 DW(2*i-1)=DP(i); DW(2*i)=DQ(i); % DV(2*i-1)=e(i); % DV(2*i)=f(i); end DW(1)=DP1; DW(2)=DV2; % DV(1)=e(1); % DV(2)=f(1); for i=1:(2*n-2) M(i)=abs(DW(i)); end if (max(M)>ac) Ja =Jacobi_netowla(G,B,e,f,n); DV=-inv(Ja)*DW; for i=1:n-1 e(i)=e(i)+DV(2*i-1); f(i)=f(i)+DV(2*i); end else break; end t=t+1; E_record=[E_record;e]; F_record=[F_record;f]; end Record=complex(E_record,F_record); V_final=Record(t,:); V_angle=angle(V_final)*180/pi; V_fuzhi=zeros(t,6); V_final_abs=abs(V_final); for j=1:t for i=1:n V_fuzhi(j,i)=sqrt(E_record(j,i).^2+F_record(j,i).^2); end end V_fuzhi_youmingzhi=V_fuzhi*220; % x=1:t; % figure(1) % plot(x,V_fuzhi(:,1),x,V_fuzhi(:,2),x,V_fuzhi(:,3),x,V_fuzhi(:,4),x,V_fuzhi(:,5),x,V_fuzhi(:,6)) % xlabel('迭代次数'); % ylabel('节点标幺值'); % legend('1','2','3','4','5','6') % %平衡节点的功率 % N=0; % for j=1:n % N=N+Y(6,j)* conj(V_final(j)); % end % S_pinghengjiedian=V_final(6)*N; % 网损计算----------------------------------------------------------------- S_ij=zeros(n); S_chuanshu=zeros(n); S_xiansun=zeros(n); I=zeros(n); % B1=[1,4,1,1,6,2,2,5; %这个矩阵是有支路两端的号数得到的； % 4,5,5,3,3,6,4,6]; % B1=[B1;R_B;X_B;B_2_B;Z]; Y1=zeros(n,n); %支路阻抗矩阵 y1=zeros(n,n); %支路导纳矩阵 for i=1:nl p=B1(1,i); q=B1(2,i); Y1(p,q)=1./B1(3,i); Y1(q,p)=Y1(p,q); y1(p,q)=B1(4,i); y1(q,p)=y1(p,q); end for i=1:n for j=1:n S_ij(i,j)=V_final(i)*V_final(i)*conj(y1(i,j))+V_final(i)*(conj(V_final(i))-conj(V_final(j)))*conj(Y1(i,j)); S_chuanshu(i,j)=V_final(i)*(conj(V_final(i))-conj(V_final(j)))*conj(Y1(i,j)); end end P1=real(S_chuanshu); Q1=imag(S_chuanshu); for i=1:n for j=1:n for k=1:nl if i==B1(1,k)&&j==B1(2,k) S_xiansun(i,j)=((P1(i,j)*P1(i,j)+Q1(i,j)*Q1(i,j))/abs(V_final(i)*V_final(i)))*B1(3,k); I(i,j)=sqrt((P1(i,j)*P1(i,j)+Q1(i,j)*Q1(i,j))/(abs(V_final(i))*abs(V_final(i)))); end end end end % 统计非零数值-------------------------------------------------------------- %S_tongji是线损的非零数值 k=1; S_tongji=[]; for i=1:n for j=1:n if S_xiansun(i,j)~=0; S_tongji=[S_tongji;k,i,j,S_xiansun(i,j)]; k=k+1; end end end S_tongji(:,4)=S_tongji(:,4)*100; % S_xiansun_shizai=zeros(1,k); % for i=1:k % S_xiansun_shizai(i)=S_tongji(i,4) %S_ij_tongji是功率传输的非零数值 S_ij_abs=abs(S_ij); k=1; S_ij_tongji=[]; for i=1:n for j=1:n if S_ij_abs(i,j)~=0; S_ij_tongji=[S_ij_tongji;k,i,j,S_ij(i,j)*100]; k=k+1; end end end %I的非零值 k=1; I_tongji=[]; for i=1:n for j=1:n if I(i,j)~=0; I_tongji=[I_tongji;k,i,j,I(i,j)]; k=k+1; end end end %平衡节点功率计算----------------------------------------------------------- sum=0; for j=1:n sum=sum+conj(Y(6,j))*conj(V_final(j)); end S_slack=V_final(6)*sum*100; % 计算线损率---------------------------------------------------------------- S_xiansun_total=0; for i=1:size(S_tongji,1) S_xiansun_total=S_xiansun_total+S_tongji(i,4); end S_total=0; for i=2:size(Ps,2) S_total=S_total-Ps(i)*100; end Xiansunlv= real(S_xiansun_total)/S_total; %各段线的线损--------------------------------------------------------------- k=1; S_ij_tongji_plus=[]; for i=1:size(S_ij_tongji,1) if real(S_ij_tongji(i,4))>0; S_ij_tongji_plus=[S_ij_tongji_plus;k,S_ij_tongji(i,2),S_ij_tongji(i,3),S_ij_tongji(i,4)]; k=k+1; end end % Lv=[]; % for i=1:size(S_ij_tongji_plus,1) % % for j=1:size(S_ij_tongji_plus,1) % if (S_ij_tongji_plus(i,2)==S_tongji(i,2))&&(S_ij_tongji_plus(i,3)==S_tongji(i,3))||(S_ij_tongji_plus(i,2)==S_tongji(i,3))&&(S_ij_tongji_plus(i,3)==S_tongji(i,2)) % Lv=[Lv;S_tongji(i,2),S_tongji(i,3),real(S_tongji(i,4))/real(S_ij_tongji_plus(i,4))]; % % elseif (S_ij_tongji_plus(i,2)==S_tongji(j,2))&&(S_ij_tongji_plus(i,3)==S_tongji(j,3))||(S_ij_tongji_plus(i,2)==S_tongji(j,3))&&(S_ij_tongji_plus(i,3)==S_tongji(j,2)) % % Lv=[Lv;S_tongji(i,2),S_tongji(i,3),real(S_tongji(j,4))/real(S_ij_tongji_plus(i,4))]; % end % end % % end %%------------------------------------------------------------------------- % fprintf('N-1检验的支路数是%d-%d\n',B1_copy(1,Ls),B1_copy(2,Ls)) % fprintf(' %.7f',V_final_abs) % fprintf(' \n') % fprintf('偏离正常率\n') % for i=1:n % fprintf(' %f%%',(V_final_abs(i)-V_fangshi1(i))*100/V_fangshi1(i)) % % end % fprintf(' \n') % for i=1:nl % fprintf('支路数是%d-%d的电流为%f,I(n-1)/I(n)=%f%%\n',I_tongji(i,2),I_tongji(i,3),I_tongji(i,4),I_tongji(i,4)*100/2.618) % end % fprintf('-------------------------------------------------------------------------\n') %%------------------------------------------------------------------------- % end

（2）雅克比矩阵求解函数

function Ja =Jacobi_netowla(G,B,e,f,n) Ja=zeros(2*(n-1),2*(n-1)); S1=0; S2=0; % PQ节点部分 for i=2:n-1 %i是PQ节点的范围 for j=1:n-1 %j是除平衡节点外的节点范围 m=2*i; t=2*j; % 非对角线的情况i!=j; if(i~=j) Ja(m-1,t-1)=-(G(i,j)*e(i)+B(i,j)*f(i)); %偏P偏e Ja(m,t)=-Ja(m-1,t-1); %偏Q偏f Ja(m-1,t)=B(i,j)*e(i)-G(i,j)*f(i); %偏P偏f Ja(m,t-1)=Ja(m-1,t); %偏Q偏e else % 对角线的情况i=j for k=1:n S1=S1+G(i,k)*e(k)-B(i,k)*f(k); S2=S2+G(i,k)*f(k)+B(i,k)*e(k); end Ja(m-1,t-1)=-S1-G(i,i)*e(i)-B(i,i)*f(i); Ja(m-1,t)=-S2+B(i,i)*e(i)-G(i,i)*f(i); Ja(m,t-1)=S2+B(i,i)*e(i)-G(i,i)*f(i); Ja(m,t)=-S1+G(i,i)*e(i)+B(i,i)*f(i); end end end %PV节点部分 S1=0; S2=0; for i=1 % i是PV节点的范围 for j=1:n-1 % j是除平衡节点外的节点 m=2*i; t=2*j; if(j~=i) Ja(m-1,t-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i)); Ja(m-1,t)=B(i,j)*e(i)-G(i,j)*f(i); Ja(m,t-1)=0; Ja(m,t)=0; else for k=1:n S1=S1+G(i,k)*e(k)-B(i,k)*f(k); S2=S2+G(i,k)*f(k)+B(i,k)*e(k); end Ja(m-1,t-1)=-S1-G(i,i)*e(i)-B(i,i)*f(i); Ja(m-1,t)=-S2+B(i,i)*e(i)-G(i,i)*f(i); Ja(m,t-1)=-2*e(i); Ja(m,t)=-2*f(i); end end end end

（3）线损计算程序

% 网损计算----------------------------------------------------------------- S_ij=zeros(n); S_chuanshu=zeros(n); S_xiansun=zeros(n); % B1=[1,4,1,1,6,2,2,5; %这个矩阵是有支路两端的号数得到的； % 4,5,5,3,3,6,4,6]; % B1=[B1;R_B;X_B;B_2_B;Z]; Y1=zeros(6,6); %支路阻抗举证 y1=zeros(6,6); %支路导纳矩阵 for i=1:nl p=B1(1,i); q=B1(2,i); Y1(p,q)=1./B1(3,i); Y1(q,p)=Y1(p,q); y1(p,q)=B1(4,i); y1(q,p)=y1(p,q); end for i=1:n for j=1:n S_ij(i,j)=V_final(i)*V_final(i)*conj(y1(i,j))+V_final(i)*(conj(V_final(i))-conj(V_final(j)))*conj(Y1(i,j)); S_chuanshu(i,j)=V_final(i)*(conj(V_final(i))-conj(V_final(j)))*conj(Y1(i,j)); end end P1=real(S_chuanshu); Q1=imag(S_chuanshu); for i=1:n for j=1:n for k=1:nl if i==B1(1,k)&j==B1(2,k) S_xiansun(i,j)=((P1(i,j)*P1(i,j)+Q1(i,j)*Q1(i,j))/abs(V_final(i)*V_final(i)))*B1(3,k); end end end end % 统计非零数值 k=1; S_tongji=[]; for i=1:n for j=1:n if S_xiansun(i,j)~=0; S_tongji=[S_tongji;i,j,k,S_xiansun(i,j)]; k=k+1; end end end

（4）三相短路的短路容量计算程序

% 短路容量计算 Ps=[7.65,4,3.5,4.2,4.2]; Qs=[0,-2.2,-1.8,-1.4,-1.5]; PQs=complex(Ps,Qs); XG=0.7; XT=0.03875; DY=zeros(6,1); DY(6)=0; for i=2:n-1 DY(i)=PQs(i); end DY(1)=3*(1/(XG+XT)); Yd=Y; for i=1:n Yd(i,i)=Y(i,i)+DY(i); end Zd=inv(Yd); Id=zeros(1,n); Sd=zeros(1,n); for f=1:n Id(f)=V_final(f)/Zd(f,f); Sd(f)=V_final(f)*conj(Id(f)); end Id_abs=abs(Id); Sd_abs=abs(Sd);