电力系统分析大作业一、设计题目本次设计题目选自课本第五章例5-8,美国西部联合电网WSCC系统的简化三机九节点系统，例题中已经给出了潮流结果，计算结果可以与之对照。

取ε=0.00001 。

二、计算步骤第一步，为了方便编程，修改节点的序号，将平衡节点放在最后。

如下图：第二步，这样得出的系统参数如下表所示：第三步，形成节点导纳矩阵。

92132 7 45683第四步，设定初值：01)0(6)0(5)0(4)0(3)0(2)0(1∠======••••••U U U U U U ；0)0(8)0(7==Q Q ，0)0(8)0(7==θθ。

第五步，计算失配功率)0(1P ∆=0，)0(2P ∆=-1.25，)0(3P ∆=-0.9，)0(4P ∆=0，)0(5P ∆=-1，)0(6P ∆=0，)0(7P ∆=1.63， )0(8P ∆=0.85；)0(1Q ∆=0.8614，)0(2Q ∆=-0.2590，)0(3Q ∆=-0.0420，)0(4Q ∆=0.6275，)0(5Q ∆=-0.1710， )0(6Q ∆=0.7101。

显然，5108614.0|},max {|-=>=∆∆εi i Q P 。

第六步，形成雅克比矩阵（阶数为14×14）第七步，解修正方程，得到：=∆)0(1θ-0.0371，=∆)0(2θ-0.0668，=∆)0(3θ-0.0628，=∆)0(4θ0.0732，=∆)0(5θ0.0191，=∆)0(6θ0.0422，=∆)0(7θ0.1726，=∆)0(8θ0.0908；=∆)0(1U 0.0334，=∆)0(2U 0.0084，=∆)0(3U 0.0223，=∆)0(4U 0.0372，=∆)0(5U 0.0266，=∆)0(6U 0.0400。

从而=)1(1θ-0.0371，=)1(2θ-0.0668，=)1(3θ-0.0628，=)1(4θ0.0732，=)1(5θ0.0191，=)1(6θ0.0422，=)1(7θ0.1726，=)1(8θ0.0908；=)1(1U 1.0334，=)1(1U 1.0084，=)1(1U 1.0223，=)1(1U 1.0372，=)1(1U 1.0266，=)1(1U 1.0400。

然后转入下一次迭代。

经三次迭代后5510101845.0|},max {|--=<⨯=∆∆εi i Q P 。

迭代过程中节点电压变化情况如下表：迭代收敛后各节点的电压和功率：最后得出迭代收敛后各支路的功率和功率损耗：三、源程序及注释由于计算流程比较简单，所以编写程序过程中没有采用模块化的形式，直接按顺序一步步进行。

disp('【节点数:】');[n1]=xlsread('input.xls','A3:A3')%节点数disp('【支路数:】');[n]=xlsread('input.xls','B3:B3')%支路数disp('【精度:】');Accuracy=xlsread('input.xls','B4:B4')%精度[branch]=xlsread('input.xls','E4:K12');[node]=xlsread('input.xls','M4:S12');Data_B1=branch;%支路参数Data_B2=node;%节点参数T1=zeros(n,2);T2=zeros(n1,3);i=sqrt(-1);format shortfor j=1:nT1(j,1)=Data_B1(j,3)+Data_B1(j,4)*1i;T1(j,2)=Data_B1(j,5)*1i;endfor j=1:n1T2(j,1)=Data_B2(j,1)+Data_B2(j,2)*1i;T2(j,2)=Data_B2(j,3)+Data_B2(j,4)*1i; endB1=zeros(n,6);B2=zeros(n1,5);for j=1:nB1(j,1)=Data_B1(j,1);B1(j,2)=Data_B1(j,2);B1(j,3)=T1(j,1);B1(j,4)=T1(j,2);B1(j,5)=Data_B1(j,6);B1(j,6)=Data_B1(j,7);endfor j=1:n1B2(j,1)=T2(j,1);B2(j,2)=T2(j,2);B2(j,3)=Data_B2(j,5);B2(j,4)=Data_B2(j,6);B2(j,5)=Data_B2(j,7);enddisp('【支路参数矩阵:】');B1 %显示支路参数矩阵disp('【节点参数矩阵:】');B2 %显示节点参数矩阵% 以上为从excel中导入初值的程序Y=zeros(n1);for i=1:nif B1(i,6)==0 %不含变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4); Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else%含有变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3);Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));endenddisp('【导纳矩阵:】');Y %显示导纳矩阵m=0;for i=1:n1if B2(i,5)==2m=m+1;endendm %PQ节点个数l=0;for i=1:n1if B2(i,5)==1l=l+1;endendl %PV节点个数Mismatch_power=zeros(l+m*2,1);for i=1:n1-1Pj=0;for j=1:n1Pj=Pj+(B2(i,3)*B2(j,3)*(real(Y(i,j))*cos(B2(i,4)-B2(j,4))+imag(Y(i,j) )*sin(B2(i,4)-B2(j,4))));endMismatch_power(i,1)=real(B2(i,1))-real(B2(i,2))-Pj;endfor k=n1:(l+m*2)Qj=0;for j=1:n1Qj=Qj+B2((k-n1+1),3)*B2(j,3)*(real(Y((k-n1+1),j))*sin(B2((k-n1+1),4)-B2(j,4))-imag(Y((k-n1+1),j))*cos(B2((k-n1+1),4)-B2(j,4)));endMismatch_power(k,1)=imag(B2((k-n1+1),1))-imag(B2((k-n1+1),2))-Qj;end% Mismatch_power %计算失配功率times=0;while(max(Mismatch_power)>Accuracy)for i=1:(n1-1)Pj=0;for j=1:n1Pj=Pj+B2(i,3)*B2(j,3)*(real(Y(i,j))*cos(B2(i,4)-B2(j,4))+imag(Y(i,j)) *sin(B2(i,4)-B2(j,4)));endMismatch_power(i,1)=real(B2(i,1))-real(B2(i,2))-Pj;endfor k=n1:(l+m*2)Qj=0;for j=1:n1Qj=Qj+B2((k-n1+1),3)*B2(j,3)*(real(Y((k-n1+1),j))*sin(B2((k-n1+1),4)-B2(j,4))-imag(Y((k-n1+1),j))*cos(B2((k-n1+1),4)-B2(j,4)));endMismatch_power(k,1)=imag(B2((k-n1+1),1))-imag(B2((k-n1+1),2))-Qj;enddisp('【当前迭代次数:】');timesdisp('【失配功率:】');Mismatch_powerJacobian=zeros(l+m*2);%雅克比矩阵7*7%————————————————————————————————————Hfor i=1:(n1-1)for j=1:(n1-1)if i==jP_H=0;for k=1:n1P_H=P_H+B2(i,3)*B2(k,3)*(real(Y(i,k))*sin(B2(i,4)-B2(k,4))-imag(Y(i,k ))*cos(B2(i,4)-B2(k,4)));endJacobian(i,i)=P_H-B2(i,3)*B2(i,3)*(0-imag(Y(i,i)));elseJacobian(i,j)=0-B2(i,3)*B2(j,3)*(real(Y(i,j))*sin(B2(i,4)-B2(j,4))-im ag(Y(i,j))*cos(B2(i,4)-B2(j,4)));endendend%————————————————————————————————————Nfor i=1:(n1-1)for j=1:mif i==jP_N=0;for k=1:n1P_N=P_N+B2(k,3)*(real(Y(i,k))*cos(B2(i,4)-B2(k,4))+imag(Y(i,k))*sin(B 2(i,4)-B2(k,4)));endJacobian(i,n1-1+i)=0-B2(i,3)*real(Y(i,i))-P_N;elseJacobian(i,n1-1+j)=0-B2(i,3)*(real(Y(i,j))*cos(B2(i,4)-B2(j,4))+imag( Y(i,j))*sin(B2(i,4)-B2(j,4)));endendend%————————————————————————————————————Kfor i=1:mfor j=1:(n1-1)if i==jP_K=0;for k=1:n1P_K=P_K+B2(i,3)*B2(k,3)*(real(Y(i,k))*cos(B2(i,4)-B2(k,4))+imag(Y(i,k ))*sin(B2(i,4)-B2(k,4)));endJacobian(n1-1+i,i)=0+B2(i,3)*B2(i,3)*real(Y(i,i))-P_K;elseJacobian(n1-1+i,j)=B2(i,3)*B2(j,3)*(real(Y(i,j))*cos(B2(i,4)-B2(j,4)) +imag(Y(i,j))*sin(B2(i,4)-B2(j,4)));endendend%————————————————————————————————————Lfor i=1:mfor j=1:mif i==jP_L=0;for k=1:n1P_L=P_L+B2(k,3)*(real(Y(i,k))*sin(B2(i,4)-B2(k,4))-imag(Y(i,k))*cos(B2(i,4)-B2(k,4)));endJacobian(n1-1+i,n1-1+i)=0-P_L+B2(i,3)*imag(Y(i,i));elseJacobian(n1-1+i,n1-1+j)=0-B2(i,3)*(real(Y(i,j))*sin(B2(i,4)-B2(j,4))-imag(Y(i,j))*cos(B2(i,4)-B2(j,4)));endendendS=zeros(l+m*2,1); %初始化电压角度变化量S=inv(Jacobian)*(0-Mismatch_power); %求解修正方程S=(Jacobian)\(0-Mismatch_power); %求解修正方程for i=1:(n1-1) %角度初值加变化量B2(i,4)=B2(i,4)+S(i,1);endfor i=1:m %电压初值加变化量B2(i,3)=B2(i,3)+S(n1-1+i,1);enddisp('【雅克比矩阵:】');Jacobian %显示雅克比矩阵% S=inv(Jacobian)times=times+1;endtimes=times-1;disp('【共计迭代次数:】');times %显示迭代次数U_It=zeros(n1,1); %初始化电压向量for i=1:n1U_It(i,1)=B2(i,3)*cos(B2(i,4))+B2(i,3)*sin(B2(i,4))*1j;endangle_It=zeros(n1,1); %将电压角度的弧度值转为角度值for i=1:n1angle_It(i,1)=B2(i,4)*180/pi;endNode_S_It=U_It.*(conj(Y)*conj(U_It)); %求解节点功率disp('【迭代收敛后各节点的电压幅值:】');Node_U_It=abs(U_It) %显示迭代收敛后各节点的电压幅值disp('【迭代收敛后各节点的电压角度:】');angle_It %显示迭代收敛后各节点的电压角度disp('【迭代收敛后各节点的功率:】');Node_S_It %显示迭代收敛后各节点的功率Branch_It=zeros(n,10);for i=1:n;if B1(i,6)==0; %不带变压器支路m=B1(i,1); %得到支路号n=B1(i,2);Branch_It(i,1)=m; %显示支路号Branch_It(i,2)=n;a=U_It(m,1)*(conj(U_It(m,1))*conj(B1(i,4))*0.5+(conj(U_It(m,1))-conj( U_It(n,1)))/conj(B1(i,3)));Branch_It(i,3)=real(a); %显示PijBranch_It(i,4)=imag(a); %显示Qijb=U_It(m,1)*B1(i,4)*0.5+(U_It(m,1)-U_It(n,1))/B1(i,3);Branch_It(i,5)=sqrt(real(b)^2+imag(b)^2); %显示Iijc=U_It(n,1)*(conj(U_It(n,1))*conj(B1(i,4))*0.5+(conj(U_It(n,1))-conj( U_It(m,1)))/conj(B1(i,3)));Branch_It(i,6)=real(c); %显示PjiBranch_It(i,7)=imag(c); %显示Qjid=U_It(n,1)*B1(i,4)*0.5+(U_It(n,1)-U_It(m,1))/B1(i,3);Branch_It(i,8)=sqrt(real(d)^2+imag(d)^2); %显示Ijie=a+c;Branch_It(i,9)=real(e); %显示线路损耗有功分量Branch_It(i,10)=imag(e); %显示线路损耗无功分量else%带变压器支路(同以上内容)m=B1(i,1);n=B1(i,2);Branch_It(i,1)=m;Branch_It(i,2)=n;a=U_It(m,1)*(conj(U_It(m,1))/conj(B1(i,3))-conj(U_It(n,1))*conj(1/(B1 (i,5)*B1(i,3))));Branch_It(i,3)=real(a);Branch_It(i,4)=imag(a);b=U_It(m,1)*(B1(i,5)-1)/B1(i,3)/B1(i,5)+(U_It(m,1)-U_It(n,1))/(B1(i,5 )*B1(i,3));Branch_It(i,5)=sqrt(real(b)^2+imag(b)^2);c=U_It(n,1)*(conj(U_It(n,1))/(conj(B1(i,5)*B1(i,5)*B1(i,3)))-conj(U_I t(m,1))*conj(1/(B1(i,5)*B1(i,3))));Branch_It(i,6)=real(c);Branch_It(i,7)=imag(c);d=U_It(n,1)*(1-B1(i,5))/B1(i,5)/B1(i,5)/B1(i,3)+(U_It(n,1)-U_It(m,1)) /B1(i,5)/B1(i,3);Branch_It(i,8)=sqrt(real(d)^2+imag(d)^2);e=a+c;Branch_It(i,9)=real(e);Branch_It(i,10)=imag(e);endenddisp('【迭代收敛后各支路的功率和功率损耗:】');Branch_It %显示迭代收敛后各支路的功率和功率损耗% %—————————————————————————————向Excel表中输出数据% Node_S_It_Real=real(Node_S_It);% Node_S_It_imag=imag(Node_S_It);% xlswrite('output.xls',Node_U_It,1,'B3');% xlswrite('output.xls',angle_It,1,'C3');% xlswrite('output.xls',Node_S_It_Real,1,'D3');% xlswrite('output.xls',Node_S_It_imag,1,'E3');% xlswrite('output.xls',Branch_It,1,'G3');程序中还有将数据从Excel表格中读入输出的xlsread和xlswrite功能，直接将数据输入到Excel表格中，可以省略将数据写在程序中或者一一输入的步骤，适用于任何节点的电力系统潮流计算。