电力系统通用潮流计算C语言程序
电力系统稳态 潮流计算程序
printf("J0[%d][%d]=%f\t",i,j,J0[i][j]);
J0[4][4]=3.6666;
//**************************************************************************************
J[i][j]=2*fi0[i];
L[i][j]=2*ei0[i];
}
else
{
H[i][j]=-ybi[i][j]*ei0[j]+ybr[i][j]*fi0[j]+Iii0[i];
N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j]+Iir0[i];
static float yd[4][4]={{0,0.01528,-0.366623,0.0192},{0.01528,0,0,0.01413},
{0.333297,0,0,0},{0.0192,0.01413,0,0}};
float ei0[4]={1,1,1,1.05};
float Max(float a[],int n)
{int i;
float max,b=0;
if(a[0]<0) {b=-a[0];}
max=b;
for(i=1;i<n;i++)
{
b=a[i];
if(a[i]<0){b=-a[i];}
if(max<b)
//*************************************************************************************
电力系统潮流计算C语言程序及说明
程序的稳定性分析
程序在不同计算机上的运行 结果是否一致。
程序运行过程中，输入数据 的变化对输出结果的影响程 度。
程序在长时间运行过程中， 输出结果是否保持稳定。
程序在处理异常情况时的表 现和稳定性。
程序的扩展性分析
代码可读性：C语言程序应具备良好的可读性，方便后续维护和修改 算法效率：C语言程序应采用高效的算法，提高计算速度 内存占用：C语言程序应合理利用内存，避免内存泄漏和不必要的内存占用 扩展性：C语言程序应具备良好的扩展性，方便添加新功能和优化性能
THANK YOU
汇报人：XX
程序的异常处理说明
异常类型：输入 错误、计算错误、 内存不足等
异常处理方式： 使用try-catch 语句进行异常捕 获和处理
异常处理流程： 当异常发生时， 程序会输出错误 信息并终止运行
异常处理结果： 确保程序在遇到 异常时能够正确 处理并给出相应 的提示信息
C语言程序应用示例
示例程序的输入数据格式
添加标题
添加标题
添加标题Βιβλιοθήκη 输入输出函数：用于数据的输入和 输出
函数：可重复使用的代码块，具有 特定的功能
C语言程序中电力系统模型的建立
定义节点和支路：根 据电力系统网络结构， 定义节点和支路，为 潮流计算做准备。
建立数学模型：根据 电力系统的物理特性 和元件参数，建立数 学模型，包括节点电 压、支路电流等。
实际运行时 间测试
程序的内存占用性能分析
内存占用情况：分 析程序运行过程中 内存的占用情况， 包括堆内存和栈内 存。
内存泄漏检测：检 查程序是否存在内 存泄漏，即程序运 行结束后未正确释 放的内存。
内存使用优化：根 据内存占用情况， 优化程序中的数据 结构或算法，降低 内存占用。
电力系统潮流计算完整c语言程序(含网损计算的最终版)
{
ia[i]=ia[i]+gY_G[n][j]*ge[j]-gY_B[n][j]*gf[j];
ib[i]=ib[i]+gY_G[n][j]*gf[j]+gY_B[n][j]*ge[j];
}
}
for(i=0,n=1;i<Bus_Num-1;i++,n++)
{
gDelta_PQ[2*i]=gDelta_P[i];
gDelta_PQ[2*i+1]=gDelta_Q[i];
}
if((fp=fopen("C:\\Documents and Settings\\Zorro\\桌面\\1\\data\\unbalance.txt","w"))==NULL)
if(gBus[n].Type==1)
gDelta_Q[i]=gDelta_Q[i]-gf[n]*(gY_G[n][j]*ge[j]-gY_B[n][j]*gf[j])+ge[n]*(gY_G[n][j]*gf[j]+gY_B[n][j]*ge[j]);
}
}
for(i=0;i<Bus_Num-1;i++)
{
gY_G[i][j]=0.0;
gY_B[i][j]=0.0;
}
for(l=0;l<Line_Num;l++)
{
i=gLine[l].No_I-1;
j=gLine[l].No_J-1;
r=gLine[l].R;
x=gLine[l].X;
电力系统潮流及短路电流计算程序
班级：姓名：学号：一、作业要求编写程序计算图1所示算例系统的潮流及三相短路电流..潮流计算：方法不限;计算系统的节点电压和相角..短路电流：4号母线发生金属性三相短路时z f=0;分别按照精确算法和近似算法计算短路电流、系统中各节点电压以及网络中各支路的电流分布;并对两种情况下的计算结果进行比较..二、电路图及参数图1 3机9节点系统表1 9节点系统支路参数表2 9节点系统发电机参数表3 9节点系统负荷参数三、计算步骤(1) 进行系统正常运行状态的潮流计算;求得(0)i U (2) 形成不含发电机和负荷的节点导纳矩阵Y N ；(3) 将发电机表示为电流源i I /i diE jx ''=和导纳i y 1/di jx '=的并联组合；节点负荷用恒阻抗的接地支路表示；形成包括所有发电机支路和负荷支路的节点导纳矩阵Y;即在Y N 中的发电机节点和负荷节点的自导纳上分别增加发电机导纳i y 和负荷导纳,LD i y *,,22LD i LDi LDiLD i i i S P jQ y V V -==； (4) 利用1Z Y-=;计算节点阻抗矩阵;从而得到阻抗矩阵中的第f 列；(5) 利用公式6-7或6-10计算短路电流；(6)利用公式6-8或6-11计算系统中各节点电压；(7)利用公式6-9计算变压器支路的电流；对输电线路利用П型等值电路计算支路电流..四、计算结果节点导纳矩阵Yn：Columns 1 through 50 -17.3611i 0 0 0 +17.3611i 00 0 -16.0000i 0 0 00 0 0 -17.0648i 0 00 +17.3611i 0 0 3.3074 -39.3089i -1.3652 +11.6041i0 0 0 -1.3652 +11.6041i 2.5528 -17.3382i0 0 0 -1.9422 +10.5107i 00 0 +16.0000i 0 0 -1.1876 + 5.9751i0 0 0 0 00 0 0 +17.0648i 0 0 Columns 6 through 90 0 0 00 0 +16.0000i 0 00 0 0 0 +17.0648i-1.9422 +10.5107i 0 0 00 -1.1876 + 5.9751i 0 03.2242 -15.8409i 0 0 -1.2820 + 5.5882i0 2.8047 -35.4456i -1.6171 +13.6980i 00 -1.6171 +13.6980i 2.7722 -23.3032i -1.1551 + 9.7843i-1.2820 + 5.5882i 0 -1.1551 + 9.7843i 2.4371 -32.1539i电压幅值：1.0400 1.0250 1.0250 1.0258 0.9956 1.0127 1.0258 1.0159 1.0324电压相角：0 0.1620 0.0814 -0.0387 -0.0696 -0.0644 0.0649 0.0127 0.0343节点有功：0.7164 1.6300 0.8500 0.0000 -1.2500 -0.9000 -0.0000 -1.0000 -0.0000节点无功：0.2705 0.0665 -0.1086 0.0000 -0.5000 -0.3000 -0.0000 -0.3500 -0.0000修正后的节点导纳矩阵Y：Columns 1 through 50 -20.6944i 0 0 0 +17.3611i 00 0 -19.3333i 0 0 00 0 0 -20.3982i 0 00 +17.3611i 0 0 3.3074 -39.3089i -1.3652 +11.6041i0 0 0 -1.3652 +11.6041i 3.8716 -17.6627i0 0 0 -1.9422 +10.5107i 00 0 +16.0000i 0 0 -1.1876 + 5.9751i0 0 0 0 00 0 0 +17.0648i 0 0 Columns 6 through 90 0 0 00 0 +16.0000i 0 00 0 0 0 +17.0648i-1.9422 +10.5107i 0 0 00 -1.1876 + 5.9751i 0 04.1321 -16.0184i 0 0 -1.2820 +5.5882i0 2.8047 -35.4456i -1.6171 +13.6980i 00 -1.6171 +13.6980i 3.7323 -23.6669i -1.1551 + 9.7843i-1.2820 + 5.5882i 0 -1.1551 + 9.7843i 2.4371 -32.1539i节点阻抗矩阵Z的第4列：0.0463 + 0.1252i0.0329 + 0.0693i0.0316 + 0.0707i0.0552 + 0.1493i0.0589 + 0.1204i0.0562 + 0.1226i0.0397 + 0.0838i0.0416 + 0.0814i0.0378 + 0.0845i精确计算结果：短路电流：模值：6.4459相角：-71.9365节点电压模值：0.1831 0.5687 0.5427 0.0000 0.1466 0.1506 0.4537 0.4463 0.4495支路电流：i j Iij1 4 0.5779-3.1264i2 7 1.3702-1.4433i3 9 0.64294-1.4808i4 5 -0.77968+1.5248i4 6 -0.6411+1.477i5 7 -0.89528+1.6436i6 9 -0.73353+1.5487i7 8 0.50734+0.10234i8 9 0.062766+0.056451i近似计算结果：短路电流：模值：6.2838相角：-69.7198节点电压模值：0.1611 0.5214 0.5157 0.0000 0.1827 0.1675 0.4227 0.4348 0.4217五、程序流程图六、程序及输入文件input_data.xls 文件：powerflow_cal.m 文件：l=9;%支路数n=9;%节点数m=6;%PQ节点数Yn=zerosn;%初始化节点导纳矩阵Y DATA1=xlsread'input_data.xls';1; %计算节点导纳矩阵Yfor k=1:li=DATA1k;1;j=DATA1k;2;R=DATA1k;3;X=DATA1k;4;B2=DATA1k;5;Yni;i=Yni;i+1i*B2+1/R+1i*X; Ynj;j=Ynj;j+1i*B2+1/R+1i*X; Yni;j=Yni;j-1/R+1i*X;Ynj;i=Ynj;i-1/R+1i*X;enddisp'节点导纳矩阵Yn：';dispYn;G=realYn;B=imagYn;DATA2=xlsread'input_data.xls';2;P=zeros1;n;Q=zeros1;n;U=ones1;n;P2:n=DATA22:n;3;Q4:n=DATA24:n;4;U1:3=DATA21:3;5;%设置节点电压初值e1=DATA21;5;e2:n=1.0;f1:n=0.0;%设置迭代次数t=0;tmax=10;while t<=tmax%计算fxa1:n=0.0;c1:n=0.0;for i=2:nfor j=1:nai=ai+Gi;j*ej-Bi;j*fj;ci=ci+Gi;j*fj+Bi;j*ej;endendfor i=2:ndeltaPi=Pi-ei*ai-fi*ci;endfor j=4:ndeltaQj=Qj-fj*aj+ej*cj;endfor k=2:3deltaU2k=Uk*Uk-ek*ek-fk*fk;endfx=deltaP2:n deltaQ4:n deltaU22:3';%计算雅克比矩阵Jfor i=2:nfor j=2:nif i~=jHi;j=-Gi;j*ei+Bi;j*fi;Ni;j=Bi;j*ei-Gi;j*fi;elseHi;j=-ai-Gi;i*ei+Bi;i*fi;Ni;j=-ci+Bi;i*ei-Gi;i*fi;endendendfor i=4:nfor j=2:nif i~=jMi;j=Bi;j*ei-Gi;j*fi;Li;j=Gi;j*ei+Bi;j*fi;elseMi;j=ci+Bi;i*ei-Gi;i*fi;Li;j=-ai+Gi;i*ei+Bi;i*fi;endendendfor i=2:3for j=2:nif i~=jRi;j=0;Si;j=0;elseRi;j=-2*ei;Si;j=-2*fi;endendendJ=H2:n;2:n N2:n;2:n;M4:n;2:n L4:n;2:n;R2:3;2:n S2:3;2:n;if maxabsfx<0.0001%输出结果break;else%求解修正方程获得dxdx=-J^-1*fx;dx=dx';e2:n=e2:n+dx1:n-1;f2:n=f2:n+dxn:2*n-1;t=t+1;endendif t>tmaxstr='潮流计算不收敛';dispstr;elsea1:n=0.0;c1:n=0.0;for i=1:nfor j=1:nai=ai+Gi;j*ej-Bi;j*fj; ci=ci+Gi;j*fj+Bi;j*ej;endendfor i=1:nUi=ei+1i*fi;ampi=absUi;argi=angleUi;Pi=ei*ai+fi*ci;Qi=fi*ai-ei*ci;enddisp'电压幅值：';dispamp;disp'电压相角：';disparg;disp'节点有功：';dispP;disp'节点无功：';dispQ;end%计算短路电流f=4;zf=0.0;%修正节点导纳矩阵Xd=DATA21:3;6;E=DATA21:3;7;for i=1:3Iii=Ei/1i*Xdi;endY=Yn;for i=1:3Yi;i=Yi;i+1/1i*Xdi;endfor j=4:nYj;j=Yj;j+-Pj+1i*Qj/Uj*Uj;enddisp'修正后的节点导纳矩阵Y：';Z=Y^-1;disp'节点阻抗矩阵Z的第4列：';dispZ:;4;%精确计算disp'精确计算结果：';U0=U;If=U0f/Zf;f+zf;amp=absIf;arg=atandimagIf/realIf;disp'短路电流：';disp'模值：';dispamp;disp'相角：';disparg;for i=1:nUi=U0i-Zi;f*If;amp=absU;enddisp'节点电压模值：';dispamp;disp'支路电流： ';str='i ''j '' Iij';dispstr;for k=1:li=DATA1k;1;j=DATA1k;2;r=DATA1k;3;x=DATA1k;4;z=r+1i*x;I=Ui-Uj/z;str=num2stri ' ' num2strj ' ' num2strI; dispstr;end%近似计算disp'近似计算结果：';U01:n=1.0;If=U0f/Zf;f+zf;amp=absIf;arg=atandimagIf/realIf;disp'短路电流：';disp'模值：';dispamp;disp'相角：';for i=1:nUi=U0i-Zi;f*If; amp=absU;enddisp'节点电压模值：'; dispamp;。
电力系统分析潮流计算C语言编程-pq分解法2
void solve(float **B,float *X,int N);/*解方程组*/
void PrtNode();/*打印输出节点参数*/
void ErrorMsg(int Flag);/*错误信息*/
int Node;/*节点数*/
int num;/*原始节点序号*/
kp=0;
for(i=0;i<NP;i++)
{
dPi=dP+i;
Yi=*(Y+i)-i;
Dltai=*(Dlta+i);
*dPi=0;
for(j=0;j<Node;j++)
{
temp=Dltai-*(Dlta+j);
if(i>j)*dPi+=*(V+j)*(Pji);
tP=*(V+j)*(Pij);
tP=*(V+i)*Yij.G-tP;
tP*=*(V+i);
tQ=*(V+j)*(Qij);
tQ-=*(V+i)*(Yij.B-Yij.B0);
tQ*=*(V+i);
}
fprintf(out,"S[%d,%d]=(%10.6f,%10.6f)\n",k+1,m+1,-tP,-tQ)
*(num+i)=k;
fscanf(in,"%d",&k);
}
if(NQ+j!=Node)ErrorMsg(4);
fprintf(out,"【节点参数表】\n");
C语言计算潮流程序
节点数：3 支路数：3 计算精度：0.00010支路1：0.0300+j0.09001┠—————□—————┨2支路2：0.0200+j0.09002┠—————□—————┨3支路3：0.0300+j0.09003┠—————□—————┨1节点1：PQ节点，S（1）=-0.5000-j0.2000节点2：PQ节点，S（2）=-0.6000-j0.2500节点3：平衡节点，U（3）=1.0000∠0.0000n=5;nl=5;isb=1;pr=0.00001;B1=[1 2 0.03i 0 1.05 0;2 3 0.08+0.3i 0.5i 1 0;2 4 0.1+0.35i 0 1 0;3 4 0.04+0.25i 0.5i 1 0;3 5 0.015i 0 1.05 1];B2=[0 0 1.05 1.05 0 1;0 3.7+1.3i 1.05 0 0 2;0 2+1i 1.05 0 0 2;0 1.6+0.8i 1.05 0 0 2;5 0 1.05 1.05 0 3];X=[1 0;2 0;3 0;4 0;5 0];na=3;Y=zeros(n);YI=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);O=zeros(1,n); for i=1:nif X(i,2)~=0;p=X(i,1);Y(p,p)=1./X(i,2);endendfor i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));YI(p,q)=YI(p,q)-1./B1(i,3);Y(q,p)=Y(p,q);YI(q,p)=YI(p,q);Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;YI(q,q)=YI(q,q)+1./B1(i,3);Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;YI(p,p)=YI(p,p)+1./B1(i,3);endG=real(Y);B=imag(YI);BI=imag(Y);for i=1:nS(i)=B2(i,1)-B2(i,2);BI(i,i)=BI(i,i)+B2(i,5);endP=real(S);Q=imag(S);for i=1:ne(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4);endfor i=1:nif B2(i,6)==2V(i)=sqrt(e(i)^2+f(i)^2);O(i)=atan(f(i)./e(i));endendfor i=2:nif i==nB(i,i)=1./B(i,i);else IC1=i+1;for j1=IC1:nB(i,j1)=B(i,j1)./B(i,i);endB(i,i)=1./B(i,i);for k=i+1:nfor j1=i+1:nB(k,j1)=B(k,j1)-B(k,i)*B(i,j1);endendendendp=0;q=0;for i=1:nif B2(i,6)==2p=p+1;k=0;for j1=1:nif B2(j1,6)==2k=k+1;A(p,k)=BI(i,j1);endendendendfor i=1:naif i==naA(i,i)=1./A(i,i);else k=i+1;for j1=k:naA(i,j1)=A(i,j1)./A(i,i);endA(i,i)=1./A(i,i);for k=i+1:nafor j1=i+1:naA(k,j1)=A(k,j1)-A(k,i)*A(i,j1);endendendendICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1;while ICT2~=0|ICT3~=0ICT2=0;ICT3=0;for i=1:nif i~=isbC(i)=0;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+BI(i,k)*sin(O(i)-O(k)));endDP1(i)=P(i)-V(i)*C(i);DP(i)=DP1(i)./V(i);DET=abs(DP1(i));if DET>=prICT2=ICT2+1;endendendNp(K)=ICT2;if ICT2~=0for i=2:nDP(i)=B(i,i)*DP(i);if i~=nIC1=i+1;for k=IC1:nDP(k)=DP(k)-B(k,i)*DP(i);endelsefor LZ=3:iL=i+3-LZ;IC4=L-1;for MZ=2:IC4I=IC4+2-MZ;DP(I)=DP(I)-B(I,L)*DP(L);endendendendfor i=2:nO(i)=O(i)-DP(i);endkq=1;L=0;for i=1:nif B2(i,6)==2C(i)=0;L=L+1;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));endDQ1(i)=Q(i)-V(i)*C(i);DQ(L)=DQ1(i)./V(i);DET=abs(DQ1(i));if DET>=prICT3=ICT3+1;endendendelse kp=0;if kq~=0;L=0;for i=1:nif B2(i,6)==2C(i)=0;L=L+1;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));endDQ1(i)=Q(i)-V(i)*C(i);DQ(L)=DQ1(i)./V(i);DET=abs(DQ1(i));endendendendNq(K)=ICT3;if ICT3~=0L=0;for i=1:naDQ(i)=A(i,i)*DQ(i);if i==nafor LZ=2:iL=i+2-LZ;IC4=L-1;for MZ=1:IC4I=IC4+1-MZ;DQ(I)=DQ(I)-A(I,L)*DQ(L);endendelseIC1=i+1;for k=IC1:naDQ(k)=DQ(k)-A(k,i)*DQ(i);endendendL=0;for i=1:nif B2(i,6)==2L=L+1;V(i)=V(i)-DQ(L);endendkp=1;K=K+1;elsekq=0;if kp~=0K=K+1;endendfor i=1:nDy(K-1,i)=V(i);endenddisp('迭代次数')disp(K);disp('每次没有达到精度要求的有功功率个数为'); disp(Np);disp('每次没有达到精度要求的无功功率个数为'); disp(Nq);for k=1:nE(k)=V(k)*cos(O(k))+V(k)*sin(O(k))*j;O(k)=O(k)*180./pi;enddisp('各节点的电压标么值E为');disp(E);disp('各节点的电压V大小');disp(V);disp('各节点的电压相角O');disp(O);for p=1:nC(p)=0;for q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));endS(p)=E(p)*C(p);enddisp('各节点的功率为');disp(S);disp('各条支路的首端功率为');for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSi(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1. /(B1(i,3)*B1(i,5))));disp(Si(p,q));enddisp('各条支路的末端功率为');for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1 ./(B1(i,3)*B1(i,5))));disp(Sj(q,p));enddisp('各条支路的功率损耗为');for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endDS(i)=Si(p,q)+Sj(q,p);disp(DS(i));endfor i=1:KCs(i)=i;for j=1:nDy(K,j)=Dy(K-1,j);endenddisp('每次迭代后各节点的电压值如图所示'); plot(Cs,Dy)xlabel('迭代次数')ylabel('电压')title('电压迭代次数曲线');题号：2节点数：4 支路数：4 计算精度：0.00010支路1：0.0200+j0.08001┠—————□—————┨3支路2：0.0400+j0.12001┠—————□—————┨4支路3：0.0500+j0.14002┠—————□—————┨4支路4：0.0400+j0.12003┠—————□—————┨4节点1：PQ节点，S(1)=-0.6000-j0.2500节点2：PQ节点，S(2)=-0.8000-j0.3500节点3：PV节点，P(3)=0.4000 V(3)=0.9500节点4：平衡节点，U(4)=1.0000∠0.0000#include "stdio.h"#include "math.h"#define M 20 /*节点数、支路数极限值*/#define N 5 /*迭代次数极限值*/int n,m,dd=0,pq=0,pv=0,tt,qy;float eps; /*节点数、支路数、对地支路数、pq节点数、pv节点数、控制打印、互联网数、精度*/struct jiedian /*节点数据结构体*/{int s; /*节点类型：1-PQ节点；2-PV节点；3-平衡节点*/ float p,q,e,f,v; /*节点的有功、无功、电压实部、虚部、电压辐值*/}jiedian[M]; /*如引用节点1的无功，则为jiedian[1].q */struct zhilu /*支路数据结构体*/{int p1,p2,s; /*支路两端节点号，支路类型：1-普通支路；2-变压器支路；3-对地支路*/float r,x,b,kt; /*支路的电阻、电抗、导纳、变压器非标准变比*/}zhilu[M];struct hulianwang /*互联网结构体*/{int num,pv; /*互联区域号,指定的pv节点号*/float p,eps2; /*规定的有功功率及其允许误差*/int count; /*每个互联区域包括的节点数*/int a[M]; /*包括的节点*/} hulianwang[M];static float G[M][M],B[M][M],G1[M][M],B1[M][M]; /*节点导纳阵矩阵*/static float ykb[2*M][2*M]; /*节点导纳阵矩阵*/static float yinzi[2*M][2*M]; /*因子表*/static float P[M][M],Q[M][M]; /*潮流计算结果*/FILE *fp1,*fp2; /*文件输入、输出指针*/void input() /*从文件in.txt 输入线路基本参数、节点数据、支路数据*/{int i,j,h; /*循环变量、节点（支路）类型*//**********************************打开电网参数文件**********************/fp1=fopen("in.txt","r");if (fp1==NULL){printf("Can not open file in.txt !\n");exit(0);}/***********读取节点数,支路数,互联网数,PQ节点数PV节点数和精度***********/fscanf(fp1,"%d,%d,%d,%f\n",&n,&m,&qy,&eps);/************************读取节点信息*************************************/for(i=1;i<=n;i++){fscanf(fp1,"%d",&h);if(h==1) /*类型h=1是PQ节点*/{pq++;fscanf(fp1,",%f,%f,%f,%f\n",&jiedian[i].p,&jiedian[i].q,&jiedian[i].e,&jiedian[i].f);jiedian[i].s=1;jiedian[i].v=sqrt(jiedian[i].e*jiedian[i].e+jiedian[i].f*jiedian[i].f);}else if(h==2) /*类型h=2是pv节点*/{pv++;fscanf(fp1,",%f,%f,%f,%f\n",&jiedian[i].p,&jiedian[i].v,&jiedian[i].e,&jiedian[i].f);jiedian[i].s=2;jiedian[i].q=0;}else /*类型h=3是平衡节点*/{jiedian[i].p=0;jiedian[i].q=0;jiedian[i].e=1;jiedian[i].f=0;jiedian[i].v=1;jiedian[i].s=3;}}/*******************************读取支路信息*************************************/for(i=1;i<=m;i++){fscanf(fp1,"%d",&h);if(h==1) /*类型h=1是普通支路*/{fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,& zhilu[i].b);zhilu[i].kt=1;zhilu[i].s=1;}if(h==2) /*类型h=2是变压器支路*/{fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,& zhilu[i].kt);zhilu[i].s=2;}if(h==3) /*类型h=3是接地支路*/{fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,& zhilu[i].b);zhilu[i].kt=1;zhilu[i].s=3;dd++;}}/*******************************读互联网信息*************************************/if(qy!=0)for(i=1;i<=qy;i++) /*输入互联网状况*/{fscanf(fp1,"%d,%d,%f,%f,%d",&hulianwang[i].num,&hulianwang[i].pv,&hulian wang[i].p,&hulianwang[i].eps2,&hulianwang[i].count);for(j=1;j<=hulianwang[i].count;j++)fscanf(fp1,",%d",&(hulianwang[i].a[j]));}fclose(fp1);/**********************************打开输出结果文件********************************/fp2=fopen("out.txt","w");if(fp2==NULL){printf("Can not open file!\n");exit(0);}fprintf(fp2,"\n**************************** 原始数据如下***********************************\n\n ");fprintf(fp2," 节点数:%2d 支路数:%2d 对地支路数:%2d PQ节点数:%2d PV节点数:%2d 精度1:%f \n",n,m,dd,pq,pv,eps);fprintf(fp2,"\n------------------------------------------------------------------------------\n\n ");for(i=1;i<=pq;i++)fprintf(fp2," 节点%2d PQ节点P[%d]=%fQ[%d]=%f\n",i,i,jiedian[i].p,i,jiedian[i].q);for(i=pq+1;i<=pq+pv;i++)fprintf(fp2," 节点%2d PV节点P[%d]=%fV[%d]=%f\n",i,i,jiedian[i].p,i,jiedian[i].v);fprintf(fp2," 节点%2d 平衡节点\n",i);fprintf(fp2,"\n-------------------------------------------------------------------------------\n\n ");for(i=1;i<=m;i++){if(zhilu[i].s==1)fprintf(fp2," 支路%2d：普通支路相关节点:%2d,%2d R=%f X=%f B=%f\n",i,zhilu[i].p1,zhilu[i].p2,zhilu[i].r,zhilu[i].x,zhilu[i].b);else if(zhilu[i].s==2)fprintf(fp2," 支路%2d：变压器支路相关节点:%2d,%2dR=%f X=%f Kt=%f\n",i,zhilu[i].p1,zhilu[i].p2,zhilu[i].r,zhilu[i].x,zhilu[i].kt);elsefprintf(fp2," 支路%2d：对地支路相关节点:%2dR=%f X=%f B=%f\n",i,zhilu[i].p1,zhilu[i].r,zhilu[i].x,zhilu[i].b);}for(i=1;i<=qy;i++){fprintf(fp2," 互联网区域:%2d 指定pv节点:%2d 规定有功功率:%8.5f 允许误差:%8.5f 包括的节点:",hulianwang[i].num,hulianwang[i].pv,hulianwang[i].p,hulianwang[i].eps2);for (j=1;j<=hulianwang[i].count;j++)fprintf(fp2,"%2d ",hulianwang[i].a[j]);fprintf(fp2,"\n");}}void youhua() /*利用节点数据、支路数据形成新编号*/ {int i,k,j,jd[N]; /*jd[N]记录节点*/struct jiedian tem; /*中间临时节点所连支路数*/for(i=1;i<=n;i++) /*对节点连接支路数赋0*/{jd[i]=0;}for(i=1;i<=n;i++) /*考虑每个节点*/{for(j=1;j<=m;j++) /*考虑每条支路*/{if((zhilu[j].p1==i||zhilu[j].p2==i)&&(zhilu[j].s!=3)) /*如果是这条支路的节点且不是接地支路*/jd[i]++;}}for(i=1;i<pq;i++) /*对pq节点按所连支路的个数进行重新编号*/{for(j=i+1;j<=pq;j++) /*从小到大排序算法*/{if(jd[i]>jd[j]){tem=jiedian[i];jiedian[i]=jiedian[j];jiedian[j]=tem;for(k=1;k<=m;k++) /*更新支路信息*/{if(zhilu[k].p1==j)zhilu[k].p1=i;else if(zhilu[k].p2==j)zhilu[k].p2=i;else if(zhilu[k].p1==i)zhilu[k].p1=j;else if(zhilu[k].p2==i)zhilu[k].p2=j;}}}}for(i=pq+1;i<pq+pv;i++) /*对pv节点按所连支路的个数进行重新编号*/{for(j=i+1;j<=pq+pv;j++) /*从小到大排序算法*/{if(jd[i]>jd[j]){tem=jiedian[i];jiedian[i]=jiedian[j];jiedian[j]=tem;for(k=1;k<=m;k++) /*更新支路信息*/{if(zhilu[k].p1==j)zhilu[k].p1=i;else if(zhilu[k].p2==j)zhilu[k].p2=i;else if(zhilu[k].p1==i)zhilu[k].p1=j;else if(zhilu[k].p2==i)zhilu[k].p2=j;}}}}fprintf(fp2,"\n**************************** 节点优化结果如下*********************************\n\n ");for(i=1;i<=n;i++){if(jiedian[i].s==1)fprintf(fp2," 节点%2d PQ节点P[%d]=%fQ[%d]=%f \n",i,i,jiedian[i].p,i,jiedian[i].q);if(jiedian[i].s==2)fprintf(fp2," 节点%2d PV节点P[%d]=%fV[%d]=%f \n",i,i,jiedian[i].p,i,jiedian[i].v);if(jiedian[i].s==3)fprintf(fp2," 节点%2d 平衡节点\n",i);}}void form_y() /*利用支路数据形成Y，注意对地支路*/ {int i,j,k;float S;for(i=1;i<=n;i++)for(j=1;j<=n;j++)G[i][j]=B[i][j]=0;for(i=1;i<=m;i++) /*节点导纳矩阵的主对角线上的导纳*/ for(j=1;j<=n;j++)if(zhilu[i].p1==j||zhilu[i].p2==j){S=zhilu[i].r*zhilu[i].r+zhilu[i].x*zhilu[i].x;if(S==0) continue;G[j][j]+=zhilu[i].r/S;B[j][j]+=-zhilu[i].x/S;if(zhilu[i].s==1) /*如果是普通支路*/{B[j][j]+=zhilu[i].b/2;}else if(zhilu[i].s==2) /*如果是普通变压器支路*/{if(zhilu[i].p1==j){G[j][j]+=(zhilu[i].r/S*(1-zhilu[i].kt))/(zhilu[i].kt*zhilu[i].kt);B[j][j]+=(-zhilu[i].x/S*(1-zhilu[i].kt))/(zhilu[i].kt*zhilu[i].kt);}else if(zhilu[i].p2==j){G[j][j]+=(zhilu[i].r/S*(zhilu[i].kt-1))/zhilu[i].kt;B[j][j]+=(-zhilu[i].x/S*(zhilu[i].kt-1))/zhilu[i].kt;}}else if(zhilu[i].s==3) /*如果是对地支路*/{B[j][j]+=zhilu[i].b;}}for(k=1;k<=m;k++) /*节点导纳矩阵非主对角线上的导纳*/ {i=zhilu[k].p1;j=zhilu[k].p2;S=zhilu[k].r*zhilu[k].r+zhilu[k].x*zhilu[k].x;if(S==0) continue;G[i][j]+=-zhilu[k].r/S;B[i][j]+=zhilu[k].x/S;if(zhilu[k].kt!=1.0){G[i][j]/=zhilu[k].kt;B[i][j]/=zhilu[k].kt;}G[j][i]=G[i][j];B[j][i]=B[i][j];}fprintf(fp2,"\n\n****************************** 节点导纳矩阵为**********************************\n");for(i=1;i<=n;i++){fprintf(fp2,"\n ");for(j=1;j<=n;j++)fprintf(fp2,"%8.5f+j%8.5f ",G[i][j],B[i][j]);}}void form_j() /*利用节点数据和Y形成J*/{float ei,fi,a=0,b=0;int i1,j1,k1,i,j,k;for(i=1;i<=2*(pq+pv)+1;i++)for(j=1;j<=2*(pq+pv)+1;j++)ykb[i][j]=0;for(i=1;i<=pq;i++)for(j=1;j<n;j++){ i1=i;j1=j;ei=jiedian[i].e;fi=jiedian[i].f;if(i!=j) /*求i!=j时的H、N、J、L*/ { ykb[2*i-1][2*j-1]=G[i1][j1]*ei+B[i1][j1]*fi; /* H */ykb[2*i-1][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi; /* N */ykb[2*i][2*j-1]=-B[i1][j1]*ei+G[i1][j1]*fi; /* J */ykb[2*i][2*j]=-G[i1][j1]*ei-B[i1][j1]*fi; /* L */}else /*求i=j时的H、N、J、K*/ { a=0;b=0;for(k=1;k<=n;k++)if(k!=i){ k1=k;a=a+G[i1][k1]*jiedian[k].e-B[i1][k1]*jiedian[k].f;b=b+G[i1][k1]*jiedian[k].f+B[i1][k1]*jiedian[k].e;}ykb[2*i-1][2*j-1]=2*G[i1][i1]*jiedian[i].e+a; /*H*/ykb[2*i][2*j]=-2*B[i1][i1]*jiedian[i].f+a; /*L*/ykb[2*i-1][2*j]= 2*G[i1][i1]*jiedian[i].f+b; /*N*/ykb[2*i][2*j-1]=-2*B[i1][i1]*jiedian[i].e-b; /*J*/}}for(i=pq+1;i<=pq+pv;i++) /*形成pv节点子阵*/for(j=1;j<n;j++){ i1=i;j1=j;ei=jiedian[i].e;fi=jiedian[i].f;if(i!=j) /*求i!=j时的H、N*/ {ykb[2*i-1][2*j-1]=G[i1][j1]*ei+B[i1][j1]*fi; /*H*/ykb[2*i-1][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi; /*N*/}else /*求i=j时的H、N、R、S*/{ a=0;b=0;for(k=1;k<=n;k++)if(k!=i){ k1=k;a+=G[i1][k1]*jiedian[k].e-B[i1][k1]*jiedian[k].f;b+=G[i1][k1]*jiedian[k].f+B[i1][k1]*jiedian[k].e;}ykb[2*i-1][2*j-1]=2*G[i1][i1]*jiedian[i].e+a; /*H*/ykb[2*i-1][2*i]=2*G[i1][i1]*jiedian[i].f+b; /*N*/ykb[2*i][2*j-1]=2*ei; /*R*/ykb[2*i][2*j]=2*fi; /*S*/}}fprintf(fp2,"\n\n\n******************************* 雅可比矩阵*************************************** \n");for(i=1;i<=2*(pq+pv);i++){fprintf(fp2,"\n ");for(j=1;j<=2*(pq+pv);j++)fprintf(fp2,"%8.5f ",ykb[i][j]);}}void form_yz() /*利用J形成因子表*/{int i,j,k;float a[2*M][2*M],L[2*M][2*M],R[2*M][2*M];float x,y,z;for(i=1;i<=2*(pq+pv);i++)for(j=1;j<=2*(pq+pv);j++){ a[i][j]=0;R[i][j]=0; L[i][j]=0;}for(i=1;i<=2*(pq+pv);i++)for(j=1;j<=2*(pq+pv);j++)a[i][j]=ykb[i][j];for(i=1;i<=2*(pq+pv);i++){for(j=1;j<=2*(pq+pv);j++){if(i>j){R[i][j]=0;x=0;for(k=1;k<=j-1;k++){x=x+L[i][k]*R[k][j];}L[i][j]=a[i][j]-x;}if(i==j){R[i][j]=1;y=0;for(k=1;k<=i-1;k++){y=y+L[i][k]*R[k][i];}L[i][i]=a[i][j]-y;}if(i<j){L[i][j]=0;z=0;for(k=1;k<=i-1;k++){z=z+L[i][k]*R[k][j];}R[i][j]=(a[i][j]-z)/L[i][i];}}}for(i=1;i<=2*(pq+pv);i++){for(j=1;j<=2*(pq+pv);j++){if(i<j)yinzi[i][j]=R[i][j];if(i==j)yinzi[i][j]=1/L[i][j];if(i>j)yinzi[i][j]=L[i][j];}}fprintf(fp2,"\n\n\n******************************** 因子表***************************************** \n");for(i=1;i<=2*(pq+pv);i++){fprintf(fp2,"\n ");for(j=1;j<=2*(pq+pv);j++)fprintf(fp2,"%8.5f ",yinzi[i][j]);}}。
电力系统潮流计算C语言程序及说明
实验目的根据所给的电力系统，编制潮流计算程序，通过计算机进行调试，最后完成一个切实可行的电力系统计算应用程序.通过自己设计电力系统计算程序使同学对电力系统分析有进一步理解，同时加强计算机实际应用能力的训练。
程序计算原理1、概述应用计算机进行电力系统计算，首先要掌握电力系统相应计算的数学模型；其次是运用合理的计算方法；第三则是选择合适的计算机语言编制计算程序。
建立电力系统计算的相关数学模型，就是建立用于描述电力系统相应计算的有关参数间的相互关系的数学方程式。
该数学模型的建立往往要突出问题的主要方面，即考虑影响问题的主要因素，而忽略一些次要因素，使数学模型既能正确地反映实际问题，又使计算不过于复杂。
运用合理的计算方法，就是要求所选用的计算方法能快速准确地得出正确结果，同时还应要求在解算过程中占用内存少，以利提高计算机的解题规模。
选择合适的语言编写程序，就是首先确定用什么计算机语言来编制程序；其次是作出计算的流程图；第三根据流程图用选择的语言编写计算程序.然后上机调试，直到语法上无错误。
本程序采用C语言进行编程。
所编制的程序难免存在逻辑错误,因此先用一个已知结果的系统作为例题进行计算.用程序计算的结果和已知结果相比较,如果结果相差甚远就要逐步分析程序的计算步骤，查出问题的出处；如果结果比较接近,则逐步分析误差来源;直到结果正确为止。
2、电力系统潮流计算的程序算法潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压(幅值及相角）、网络中的功率分布及功率损耗等.目前计算机潮流计算的方法主要有牛顿—拉夫逊算法和PQ分解法。
牛顿-拉夫逊算法是数学上求解非线形方程组的有效方法,具有较好的收敛性,曾经是潮流计算中应用比较普遍的方法。
PQ快速分解法是从牛顿—拉夫逊算法演变而来的，是将纯数学的牛顿—拉夫逊算法与电力系统具体特点相结合并进行简化与改进而得出的。
PQ快速分解法比牛顿—拉夫逊算法大大提高了计算速度和节省了内存，故而本程序以PQ快速分解法进行潮流计算.1）形成节点导纳矩阵（1)自导纳的形成对节点i其自导纳Y ii是节点i以外的所有节点都接地时节点i对地的总导纳。
电力系统潮流计算的C语言实现
//////////////////////////////////////////////////////////////////////// PQ分解法潮流////文件输入格式：节点总数n(包括联络节点),支路数zls ////节点数(发电机和负荷)nb,接地电抗数mdk,迭代精度eps // //考虑负荷静特性标志kk2(0考虑),平衡节点号，优化标志(0不优化) ////最大迭代次数it1,支路左右节点号izl[],jzl[],支路电阻zr[],电抗zx[] ////支路容纳zyk[],节点号nob[]及标志nobt[](0-PQ,-1-PV) ////发电机和负荷有功、无功pg[],qg[],pl[],ql[] ////电压v0[](pv节点输入实际值，PQ节点任输入一值) // //电抗节点号idk[],电抗值dkk[] ////////////////////////////////////////////////////////////////////////#include "math.h"#include "stdio.h"#define NS 2000 //最大节点数#define NS2 NS * 2#define NS4 1000 //NS4、NS必须大于2*zls。
#define ZS 3000 //最大支路数#define ZS2 ZS * 2#define DKS 200 //最大电抗器数#define N2 ZS * 4#define N3 ZS * 8 + NS * 4FILE *fp1, *fp2;char inname[12], outname[12];// fp1输入数据文件指针fp2输出文件指针// inname[]输入数据文件名outname[]输出数据文件名int n, zls, nb, mdk, mpj, bnsopton, it1, dsd, kk2, nzls;// 节点总数n(包括联络节点) 支路数(回路数)zls 节点数nb(发电机和负荷) // 接地电抗数mdk 精度eps 平衡节点号mpj// 节点优化(标志)bnsopton(=0节点不优化，!=0节点优化)// 最大迭代次数it1 最低电压或最大功率误差节点号dsd// 负荷静特性标志(=0考虑负荷静特性)// 支路数(双回线算一条支路)int izl[ZS], jzl[ZS], idk[DKS], yds[NS], ydz[NS], iy[ZS2];// izl[],jzl[],idk[]:分别存放左、右节点号和电抗器节点号。
C语言进行潮流计算
C语言进行潮流计算C语言是一种通用的编程语言，被广泛用于开发各种应用程序。
虽然它的历史可以追溯到上个世纪70年代，但C语言在当今仍然非常流行。
它被广泛用于系统编程、嵌入式系统开发、游戏开发和科学计算等领域。
本文将重点介绍C语言的潮流计算。
潮流计算是一种用于分析电力系统稳态和动态行为的工具和技术。
它可以用于预测电力系统中潮流（电压和电流）的分布和变化。
潮流计算的目标是找到系统中各个节点的电压和相应的功率流向。
这对于电力系统的规划和运行至关重要。
C语言可以通过编写相应的程序来实现潮流计算。
在这个程序中，需要定义一组节点和支路，以及相应的导纳矩阵（节点之间的电导和电纳）。
导纳矩阵描述了电力系统中节点之间的电压和电流关系。
一旦定义了这些参数，就可以使用C语言编写算法来解决潮流计算问题。
潮流计算通常使用迭代的方法来求解。
在每一次迭代中，系统的节点电压和功率分布都会进行更新，直到达到稳定状态。
C语言提供了丰富的控制结构和数据类型，使得编写这样的迭代算法变得相对容易。
例如，可以使用for循环来实现迭代过程，同时使用浮点数类型来存储节点电压和功率分布。
潮流计算的核心是求解节点电压的方程组。
这个方程组通常是一个非线性方程组，需要使用牛顿-拉夫逊方法或高斯-赛德尔方法等迭代算法进行求解。
C语言提供了一些数值计算库，如GSL（GNU科学库），可以用于实现这些迭代算法。
这些库提供了各种数值计算函数，如求解非线性方程组、矩阵运算和插值等，可以大大简化潮流计算程序的编写。
除了求解节点电压方程组，潮流计算还需要考虑各种系统参数和条件。
例如，支路的参数（如电阻和电感）、发电机的容量和调节策略、负荷的需求和变化等。
C语言可以通过定义和存储这些参数，并在计算过程中引用它们来实现对这些条件的考虑。
此外，潮流计算还需要处理一些特殊情况，如无源系统和有环路系统。
无源系统是指没有外部电压源的电力系统，只由发电机和负荷组成。
有环路系统是指电力系统中存在环路，这会导致方程组无解。
Matlab及C语言在潮流计算运用
Matlab 及C 语言在潮流计算中的运用陈洲<三峡大学电气信息学院20071096班）摘要：潮流计算是研究电力系统稳态运行情况的一种基本电气计算，常规潮流计算的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压<幅值及相角）、网络中的功率分布以及功率损耗等。
潮流计算的结果是电力系统稳定计算和故障分析的基础。
关键词：电力系统分析 潮流计算 牛顿-拉夫逊法 C 语言 Matlab 一，潮流计算算法原理：牛顿—拉夫逊法的基本原理牛顿-拉夫逊法是一种求解非线性方程的数值解法，因为便于编写程序用计算机求解，应用较广。
下面以一元非线性代数方程的求解为例，来说明牛顿-拉夫逊法的基本思想。
设欲求解的非线性代数方程为 f(x>=o 设方程的真实解为x*，则必有f(x*>=0。
用牛顿-拉夫逊法求方程真实解x*的步骤如下： 首先选取余割合适的初始估值x°作为方程f(x>=0的解，若恰巧有f(x°>=0，则方程的真实解即为x*= x°若f(x°>≠0，则做下一步。
取x¹=x°+Δx°为第一次的修正估值，则 f(x¹>=f(x°+Δx°> 其中Δx°为初始估值的增量，即Δx°=x¹-x°。
设函数f(x>具有任意阶导数，即可将上式在x°的邻域展开为泰勒级数，即： f(x¹>=f(x°+Δx°>=f(x°>+f'(x°>Δx°+[f''(x°>(Δx°>2]/2+…若所取的|Δx°|足够小，则含(Δx°>²的项及其余的一切高阶项均可略去，并使其等于零，即： f(x¹>≈f(x°>+f'(x°>Δx°=0 故得 Δx°=-f(x°>/f'(x°> 从而 x¹= x°-f(x°>/f'(x°> 可见，只要f'(x°>≠0，即可根据上式求出第一次的修正估值x¹，若恰巧有f(x¹>=0，则方程的真实解即为x*=x¹。
电力系统通用潮流计算C语言程序
#include <iostream>#include <fstream>#include<iomanip>#include<math.h>using namespace std;//节点号类型负荷有功负荷无功母线数据（类型1＝PV节点，2＝PQ节点，3＝平衡节点）struct BUS{int busno;int type;float Pd;float Qd;};//发电机数据节点号有功发电电压幅值struct Generator{int busno;float Pg;float Vg;};//支路信息节点I 节点J R X B/2 kstruct Line{int busi;int busj;float R;float X;float B;float k;};//deltaP deltaQ deltaV^2//void fun1(double YG[][50],double YB[][50],double e[],double f[],int type[],int N,double W[],double P[],double Q[],double V[]){double dP=0,dQ=0,dV=0;int i,j;for(i=0;i<N-1;i++){double A=0,B=0;for(j=0;j<N;j++){A+=YG[i][j]*e[j]-YB[i][j]*f[j];B+=YG[i][j]*f[j]+YB[i][j]*e[j];}dV=V[i]*V[i]-e[i]*e[i]-f[i]*f[i];dP=P[i]-e[i]*A-f[i]*B;W[2*i]=dP;dQ=Q[i]-f[i]*A+e[i]*B;if(type[i]==1)W[2*i+1]=dQ;else W[2*i+1]=dV;}}//Jacobi矩阵//void Jacobi(double YG[][50],double YB[][50],double e[50],double f[50],int type[50],int N ,double Ja[100][101]){int i,j;for(i=0;i<N;i++){for(j=0;j<N;j++){if(i!=j){if(type[i]==1){Ja[2*i][2*j]=-(YG[i][j]*e[i]+YB[i][j]*f[i]);Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i];Ja[2*i+1][2*j]=Ja[2*i][2*j+1];Ja[2*i+1][2*j+1]=-Ja[2*i][2*j];}else {Ja[2*i][2*j]=-YG[i][j]*e[i]+YB[i][j]*f[i];Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i];Ja[2*i+1][2*j+1]=Ja[2*i+1][2*j]=0;}}else {double a[50]={0},b[50]={0};for(int k=0;k<N;k++){a[i]+=(YG[i][k]*e[k]-YB[i][k]*f[k]);b[i]+=(YG[i][k]*f[k]+YB[i][k]*e[k]);Ja[2*i][2*j]=-a[i]-YG[i][i]*e[i]-YB[i][i]*f[i];Ja[2*i][2*j+1]=-b[i]+YB[i][i]*e[i]-YG[i][i]*f[i];if(type[i]==1){Ja[2*i+1][2*j]=b[i]+YB[i][i]*e[i]-YG[i][i]*f[i];Ja[2*i+1][2*j+1]=-a[i]+YG[i][i]*e[i]+YB[i][i]*f[i];}else {Ja[2*i+1][2*j]=-2*e[i];Ja[2*i+1][2*j+1]=-2*f[i];}}}}}}//高斯消元法解方程组函数//void gauss(double a[][101],int n){int i,j,k;double c;for(k=0;k<n-1;k++) {c=a[k][k];for(j=k;j<=n;j++) a[k][j]/=c;for(i=k+1;i<n;i++) {c=a[i][k];for(j=k;j<=n;j++) a[i][j]-=c*a[k][j];}}a[n-1][n]/=a[n-1][n-1];for(k=n-2;k>=0;k--)for(j=k+1;j<n;j++) a[k][n]-=a[k][j]*a[j][n];}void main(){ifstream fin;int N=0,GS=0,LD=0,ZLs=0; //节点数发电机数负荷数支路数// BUS *B;Generator *G;Line *L;//从文本中读入原始数据到数组中//fin.open("C:\\data.txt");if(!fin){cout<<"输入数据文件不存在！"<<endl;getchar();}int m1[50]={0},m2[50]={0};float m3[50],m4[50],m5[50],m6[50];int i,j,l;for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m2[i]>>m3[i]>>m4[i];N++;}B =new BUS[N];for (i=0;i<N;i++){B[i].busno=m1[i];B[i].type=m2[i];B[i].Pd=m3[i];B[i].Qd=m4[i];}for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m4[i]>>m3[i];GS++;}G =new Generator[GS];for (i=0;i<GS;i++){G[i].busno=m1[i];G[i].Pg=m4[i];G[i].Vg=m3[i];}for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m2[i]>>m3[i]>>m4[i]>>m5[i]>>m6[i];ZLs++;}L =new Line[ZLs];for (i=0;i<ZLs;i++){L[i].busi=m1[i];L[i].busj=m2[i];L[i].R=m3[i];L[i].X=m4[i];L[i].B=m5[i];L[i].k=m6[i];}LD=N-GS;fin.close();//节点导纳矩阵形成//double YB[50][50],YG[50][50],BB[50][50],K[50][50];for(i=0;i<N;i++){for(j=0;j<N;j++){YB[i][j]=0;YG[i][j]=0;BB[i][j]=0;K[i][j]=1;}}for (l=0;l<ZLs;l++){i=L[l].busi-1;j=L[l].busj-1;K[i][j]=L[l].k;BB[i][j]=BB[j][i]=L[l].B;YG[i][j]=YG[j][i]=L[l].R/(L[l].R*L[l].R+L[l].X*L[l].X);YB[i][j]=YB[j][i]=-L[l].X/(L[l].R*L[l].R+L[l].X*L[l].X);}for(i=0;i<N;i++){for(j=i;j<N;j++) {K[i][j]=K[j][i];K[j][i]=1;}for(j=0;j<N;j++){if(i!=j){YG[i][i]=YG[i][i]+(YG[i][j]*K[i][j]*K[i][j]);YB[i][i]=YB[i][i]+(YB[i][j]*K[i][j]*K[i][j]+BB[i][j]);}}}//修正后//for (l=0;l<ZLs;l++){i=L[l].busi-1;j=L[l].busj-1;K[i][j]=L[l].k;YG[i][j]=-YG[i][j]*K[i][j];YG[j][i]=YG[i][j];YB[i][j]=-YB[i][j]*K[i][j];YB[j][i]=YB[i][j];}int type[50]={0};for(i=0;i<N;i++){type[i]=B[i].type;}//PQV的获得//double P[50],Q[50],V[50];for(i=0;i<N;i++){P[i]=0;Q[i]=0;V[i]=0;P[i]=-B[i].Pd;Q[i]=-B[i].Qd;}for (i=0;i<GS;i++){P[G[i].busno-1]=G[i].Pg;V[G[i].busno-1]=G[i].Vg;}// 求A=e+f//double e[50]={0},f[50]={0};double C[100]={0},D[100]={0};for(i=0;i<N;i++){if(V[i]==0){C[2*i]=1;}else C[2*i]=V[i];}double W[100]={0},Ja[100][101]={0};//调用Jacobi函数和高斯函数//for(int t=1;t<10;t++){for(i=0;i<2*N-2;i++){e[i]=C[2*i];f[i]=C[2*i+1];}fun1(YG,YB,e,f,type,N,W,P,Q,V);double it=fabs(W[0]);for(i=1;i<2*N-2;i++){if (it<fabs(W[i])) {it=fabs(W[i]);j=i;}}//中间迭代过程//cout<<setw(10)<<"迭代次数"<<setw(20)<<"最大的功率误差"<<setw(8)<<"节点号"<<endl;cout<<setw(10)<<t<<setw(20)<<it<<setw(8)<<j/2+1<<endl;if (it<0.00001) break;Jacobi(YG,YB,e,f,type,N,Ja);for(i=0;i<2*N-2;i++){Ja[i][2*N-2]=W[i];}//高斯消元法解方程//gauss(Ja,2*N-2);for(i=0;i<2*N-2;i++){D[i]=-Ja[i][2*(N-1)];C[i]+=D[i];}}//平衡节点//for(i=0;i<N;i++){double a=0,b=0;for(int j=0;j<N;j++){a+=(YG[i][j]*e[j]-YB[i][j]*f[j]);b+=(YB[i][j]*e[j]+YG[i][j]*f[j]);}P[i]=e[i]*a+f[i]*b;Q[i]=f[i]*a-e[i]*b;}//支路//double PZL[100][101]={0},QZL[100][101]={0},pr[100][101]={0},qx[100][101]={0}; double x1=0,x2=0,y1=0,y2=0,I2=0;for(int k=0;k<ZLs;k++){i=L[k].busi-1;j=L[k].busj-1;x1=e[i]/L[k].k-e[j];y1=f[i]/L[k].k-f[j];x2=-e[i]*YG[i][j]-f[i]*YB[i][j];y2=-f[i]*YG[i][j]+e[i]*YB[i][j];QZL[i][j]=(x1*y2-x2*y1);PZL[i][j]=(x1*x2+y1*y2);I2=(PZL[i][j]*PZL[i][j]+QZL[i][j]*QZL[i][j])/(e[i]*e[i]+f[i]*f[i]);pr[i][j]=I2*L[k].R;qx[i][j]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[i][j]+=(e[i]*e[i]+f[i]*f[i])*(-L[k].B);x1=e[j]*L[k].k-e[i];y1=f[j]*L[k].k-f[i];x2=-e[j]*YG[j][i]-f[j]*YB[j][i];y2=-f[j]*YG[j][i]+e[j]*YB[j][i];QZL[j][i]=(x1*y2-x2*y1);PZL[j][i]=(x1*x2+y1*y2);I2=(PZL[j][i]*PZL[j][i]+QZL[j][i]*QZL[j][i])/(e[j]*e[j]+f[j]*f[j]);pr[j][i]=I2*L[k].R;qx[j][i]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[j][i]+=(e[j]*e[j]+f[j]*f[j])*(-L[k].B);}//全网数据//int high=1,low=1;double PG=0,PL=0,Prr=0,Vh=sqrt(e[0]*e[0]+f[0]*f[0]),Vl=sqrt(e[0]*e[0]+f[0]*f[0]); for(k=0;k<N;k++){Prr+=P[k];if(B[k].type==1) PL+=B[k].Pd;else PG+=P[k];if(sqrt(e[k]*e[k]+f[k]*f[k])>Vh){Vh=sqrt(e[k]*e[k]+f[k]*f[k]);high=k+1;}if(sqrt(e[k]*e[k]+f[k]*f[k])<Vl){Vl=sqrt(e[k]*e[k]+f[k]*f[k]);low=k+1;}}//输出数据到文件databak.txt//ofstream fout;fout.open("C:\\databak.txt");fout<<"节点"<<endl;fout<<setw(8)<<"节点号"<<setw(16)<<"V"<<setw(16)<<"弧度"<<setw(16)<<"发电P"<<setw(16)<<"发电Q"<<setw(16)<<"负荷P"<<setw(16)<<"负荷Q"<<endl;for(i=0;i<LD;i++){fout<<setw(8)<<i+1<<setw(16)<<sqrt(e[i]*e[i]+f[i]*f[i])<<setw(16)<<atan2(f[i],e[i])*180/3 .14159<<setw(16)<<0<<setw(16)<<0<<setw(16)<<B[i].Pd<<setw(16)<<B[i].Qd<<endl;}for(j=0;j<GS;j++){i=G[j].busno-1;fout<<setw(8)<<i+1<<setw(16)<<V[i]<<setw(16)<<atan2(f[i],e[i])*180/3.14159<<setw(16)<<P[ i]<<setw(16)<<Q[i]<<setw(16)<<0<<setw(16)<<0<<endl;}fout<<"支路"<<endl;fout<<setw(4)<<"i"<<setw(4)<<"j"<<setw(10)<<"i_j有功"<<setw(10)<<"i_j无功"<<setw(10)<<"j_i有功"<<setw(12)<<"j_i无功"<<setw(12)<<"有功损耗"<<setw(12)<<"无功损耗"<<endl;for (k=0;k<ZLs;k++){i=L[k].busi-1;j=L[k].busj-1;fout<<setw(4)<<L[k].busi<<setw(4)<< L[k].busj <<setw(10)<<PZL[i][j] <<setw(10)<<QZL[i][j]<<setw(10)<< PZL[j][i] <<setw(12)<<QZL[j][i]<<setw(12)<<pr[i][j]<<setw(12)<<qx[i][j]<<endl;}fout<<"全网数据"<<endl;fout<<setw(14)<<"发电有功"<<setw(14)<<"负荷有功"<<setw(14)<<"有功损耗"<<setw(14)<<"最高电压"<<setw(14)<<"节点号"<<setw(14)<<"最低电压"<<setw(14)<<"节点号"<<endl;fout<<setw(14)<<PG<<setw(14)<<PL<<setw(14)<<Prr<<setw(14)<<Vh<<setw(14)<<high<<set w(14)<<Vl<<setw(14)<<low<<endl;fout.close();}。
C语言进行潮流计算
电力系统课程设计C语言潮流计算学院：电气工程班级：电092 班学号：0912002020 学生姓名：闵凯2013.3.7电力系统的潮流计算是对电力系统分析的最基本步骤也是最重要的步骤，是指在一定的系统结构和运行条件下，确定系统运行状态的计算，也即是对各母线（节点）电压，各元件（支路）传输电线或功率的计算。
通过计算出的节点电压和功率分布用以检查系统各元件是否过负荷，各点电压是否合理，以及功率损耗等。
即使对于一个简单的电力系统，潮流计算也不是一件简单就可以完成的事，其运算量很大，因此如果对于一个大的、复杂的电网来说的话，由于其节点多，分支杂，其计算量可想而知，人工对其计算也更是难上加难了。
特别是在现实生活中，遇到一个电力系统不会像我们期望的那样可以知道它的首端电压和首端功率或者是末端电压和末端功率，而是只知道它的首端电压和末端功率，更是使计算变的头疼万分。
为了使计算变的简单，我们就可以利用计算机，用C语言编程来实现牛顿-拉夫逊（Newton-Raphson）迭代法，最终实现对电力系统潮流的计算。
用牛顿-拉夫逊迭代法进行电力系统潮流计算的相关概念1.节点导纳矩阵如图所示的电力网络，将节点i和j的电压用U i和U j表示,它们之间的支路导纳表示为y j,那么有基尔霍夫电流定律可知注入接点I的电流I i（设流入节点的电流为正）等于离开节点I的电流之和，因此有y u上式也可以写为n-zj =0y j u=Y HI i:I = YU-y j =Yj则可将（1-2）改写为：Y j U j I=1,2,…,n.(1-1)(1-2)(1-3)(1-4)' y iJ(U i-u J)I i 八I ij其中Y 为节点导纳矩阵，也称为稀疏的对称矩阵，它是 nM 阶方阵。
对角元 Yu 称为自 导纳，它等于与该节点I 直接相连的所有支路导纳总和；非对角元 Y j ( i 并)称为互导纳或 转移导纳，它等于连结节点 I , j 支路导纳的负数，且有 Y j =Y ji ，当节点I , j 之间没有支路 直接相连时，Y j =Y ji =O 。
电力系统潮流运算
/*此程序为用牛顿拉夫逊法的直角坐标形式求电力线路的潮流计算程序*//*此程序为一通用程序,按给定的题目输入数据到input文件中，就可求出相关节点的电压相量，支路功率，网络损耗，输电效率平衡节点的功率，pv节点的无功*//*考虑了线路对地导纳支路、平行支路、非标准变比变压器支路、还有无功越限的处理等4种较为复杂的情况*//*本程序通用性很强*/#include <stdio.h>#include<math.h>#include<stdlib.h>void main(){int bus,line,pQ,pv,blzl,transform;/*分别表示节点数，支路数，pQ节点数，pv节点数，对地支路数，变压器数*/int i,j,k,m,n,circle,num,w,r;/*循环变量*/int type[100],node[100],head[100],end[100],branch[100];/*分别表示:节点类型，节点号，支路的首节点，支路的尾节点，支路号*/int t1[100],t2[100];/*中间数组*/float precise;/*计算精度*/float real[100][100],image[100][100],R[100],X[100];/*节点导纳的实部，节点导纳的虚部，支路上的电阻，支路上的电抗*/float p[100],Q[100],p0[100],Q0[100];/*节点的有功功率，节点的无功功率，各节点的给定有功功率初值，各节点的给定无功功率初值*/float e[100],f[100],e0[100],f0[100];/*分别为各节点电压的实部，虚部，给定电压的初始实部和虚部*/float ub[200],ubp[100],ubQ[100],ubu2[100],a[100],b[100];/*分别表示不平衡量，有功不平衡量，无功不平衡量，电压不平衡量,a[]和b[]是中间数组*/floatH[100][100],N[100][100],B[100][100],J[100][100],L[100][100],R1[100][100],S[100][100];/*雅克比矩阵各元素变量*/float a1[100][100],b1[100][100],du[100];/*中间数组和电压修正量*/floatt,s,E[100][100],Ln[100][100],Un[100][100],Bn[100][100],zgjz[100][200],Eu[100][100];/*求逆矩阵的过程中的各数组*/float dumax,ps,Qs,pline[100][100],Qline[100][100],dpline[100][100],dQline[100][100];/*最大电压修正量，平衡节点的有功和无功，支路有功和无功，支路损耗的有功和无功*/ float dpsum,dQsum,efficiency,p1,p2;/*网络总的有功损耗，网络总的无功损耗，输电效率，网络中各节点注入的正的有功的和，网络中各节点注入负的有功的绝对值和*/ float duidi[100],ratio[100],tt1,tt2;/*对地支路的导纳和变压器变比，还有中间变量tt1和tt2*/int x,y,d[100],flag[100];/*flag用来判断平行支路数*/FILE *fp1,*fp2;/*指针变量*/float ee[100],ff[100];/*中间变量*/float Qmax;/*Qmax为无功上限*//*以上为对程序中用到的变量的定义*/fp1=fopen("input.txt","r");/*打开input.txt文件读入数据*/fp2=fopen("output.txt","w");/*打开out.txt文件写入数据*/fprintf(fp2,"\n 华北电力大学\n\n");fprintf(fp2," 电力系统潮流上机班级:电气化0805 姓名:张同学号:200801000528\n");fprintf(fp2,"\n ***************原始数据***************** \n");fprintf(fp2,"============================================================= ======================\n\n");fscanf(fp1,"%d%d%d%d%d%d%f%f",&bus,&line,&pQ,&pv,&blzl,&transform,&precise,&Qmax );fprintf(fp2," 节点数:%d 支路数:%d pQ节点数:%d pv节点数:%d\n 并联对地之路数:%d 变压器数:%d 要求精度:%f \n",bus,line,pQ,pv,blzl,transform,precise);fprintf(fp2,"------------------------------------------------------------------------------------\n");/*下边开始从input.txt文件读入数据*/for(i=1;i<=bus;i++){fscanf(fp1,"%d%d",&node[i],&type[i]);if(type[i]==1)/*type[i]=1,说明第i节点是pQ节点*/fscanf(fp1,"%f%f%f%f",&p0[i],&Q0[i],&e[i],&f[i]);else if(type[i]==2)/*type[i]=1,说明第i节点是pv节点*/fscanf(fp1,"%f%f%f%f%f",&p0[i],&e[i],&f[i],&e0[i],&f0[i]);else if(type[i]==3)/*type[i]=1,说明第i节点是平衡节点*/fscanf(fp1,"%f%f",&e[i],&f[i]);}for(i=1;i<=bus;i++)/*bus表示节点数*/{if(type[i]==1)fprintf(fp2," 节点%d pQ节点%f %f\n",node[i],p0[i],Q0[i]);else if(type[i]==2)fprintf(fp2," 节点%d pv节点%f %f\n",node[i],p0[i],e0[i]);elsefprintf(fp2," 节点%d 平衡节点%f %f\n",node[i],e[i],f[i]);}fprintf(fp2,"------------------------------------------------------------------------------------\n");for(i=1;i<=line;i++){fscanf(fp1,"%d%d%d%f%f%f%d%f",&branch[i],&head[i],&end[i],&R[i],&X[i],&duidi[i],&flag[ i],&ratio[i]);/*读支路信息*/fprintf(fp2," 支路%d %d %d R[%d]=%fX[%d]=%f\n",branch[i],head[i],end[i],branch[i],R[i],branch[i],X[i]); /*打印出给定题目的支路信息*/}/*head里存放支路的首节点，end里存放支路的尾节点*/fprintf(fp2,"============================================================= ======================\n\n");fprintf(fp2,"\n *********************潮流计算开始***************** \n");fclose(fp1);/*关闭文件input.txt文件*//*下边求节点导纳阵开始*//*下面是考虑平行之路时的支路阻抗计算*/for(i=1;i<=line;i++){if(flag[i]==2){R[i]=R[i]/2;X[i]=X[i]/2;duidi[i]=duidi[i]*2;}}/*考虑平行之路的阻抗结束*/for(i=1;i<=bus;i++)for(j=1;j<=bus;j++){real[i][j]=0;/*real是节点导纳的实部*/image[i][j]=0;/*image是节点导纳的虚部*/}for(i=1;i<=line;i++)/*line是总的支路数*/{real[(head[i])][(end[i])]=(-1)*R[i]/(R[i]*R[i]+X[i]*X[i]);image[(head[i])][(end[i])]=X[i]/(R[i]*R[i]+X[i]*X[i]);real[(end[i])][(head[i])]=real[(head[i])][(end[i])];image[(end[i])][(head[i])]=image[(head[i])][(end[i])];for(j=1;j<=bus;j++)if((j==head[i])||(j==end[i])){real[j][j]+=(-1)*real[(head[i])][(end[i])];image[j][j]+=(-1)*image[(head[i])][(end[i])];}}/*下边考虑对地之路*/for(i=1;i<=line;i++){image[head[i]][head[i]]+=duidi[i];image[end[i]][end[i]]+=duidi[i];}/*对地之路考虑结束*//*下边是考虑变压器的非标准变比的程序*/for(i=1;i<=line;i++){tt1=real[(head[i])][(end[i])];tt2=image[(head[i])][(end[i])];real[(head[i])][(end[i])]=real[(head[i])][(end[i])]/fabs(ratio[i]);image[(head[i])][(end[i])]=image[(head[i])][(end[i])]/fabs(ratio[i]);real[(end[i])][(head[i])]=real[(head[i])][(end[i])];image[(end[i])][(head[i])]=image[(head[i])][(end[i])];if(ratio[i]>0){real[head[i]][head[i]]=real[head[i]][head[i]]-tt1*(1-ratio[i])/(ratio[i]*ratio[i])-real[(head[i])][(end[i ])]+tt1;image[head[i]][head[i]]=image[head[i]][head[i]]-tt2*(1-ratio[i])/(ratio[i]*ratio[i])-image[(head[i])] [(end[i])]+tt2;}if(ratio[i]<0){real[end[i]][end[i]]=real[(end[i])][(end[i])]-tt1*(1+ratio[i])/(ratio[i]*ratio[i])-real[(head[i])][(end[i] )]+tt1;image[end[i]][end[i]]=image[(end[i])][(end[i])]-tt2*(1+ratio[i])/(ratio[i]*ratio[i])-image[(head[i])] [(end[i])]+tt2;}}/*考虑变压器的非标准变比的程序结束*/fprintf(fp2,"========================================================= =================================\n\n");fprintf(fp2," 节点导纳阵为:Y=\n");/*打印节点导纳阵*/for(i=1;i<=bus;i++){for(j=1;j<=bus;j++)fprintf(fp2," %f+j%f",real[i][j],image[i][j]);fprintf(fp2,"\n");}/*节点导纳阵已经求出，并已经打印*/fprintf(fp2,"========================================================= ==================================\n\n");for (i=1;i<=bus;i++){ ee[i]=e[i];ff[i]=f[i];}loop : for (i=1;i<=bus;i++){ e[i]=ee[i];f[i]=ff[i];}fprintf(fp2,"=====================================迭代过程==============================================\n\n");fprintf(fp2,"迭代过程中的电压和功率初值为:\n");for(i=1;i<=bus;i++){if(type[i]==1)fprintf(fp2," 节点%d e=%f f=%f p=%f Q=%f\n",node[i],e[i],f[i],p0[i],Q0[i]);else if(type[i]==2)fprintf(fp2," 节点%d e=%f f=%f p=%f\n",node[i],e[i],f[i],p0[i]);}fprintf(fp2,"------------------------------------------------------------------------------------\n");circle=1;/*circle为迭代的次数*//*下边迭代开始*/do{fprintf(fp2,"第%d 次迭代\n",circle);for(i=1;i<=bus;i++){if(type[i]==1){p[i]=0;Q[i]=0;}else if(type[i]==2)p[i]=0;}/*下边是求修正方程的不平衡列向量*/fprintf(fp2,"\n修正方程的不平衡量为:\n");m=1;n=1;for(i=1;i<=bus;i++){if(type[i]==1){for(j=1;j<=bus;j++){p[i]+=(e[i]*(real[i][j]*e[j]-image[i][j]*f[j])+f[i]*(real[i][j]*f[j]+image[i][j]*e[j]));Q[i]+=(f[i]*(real[i][j]*e[j]-image[i][j]*f[j])-e[i]*(real[i][j]*f[j]+image[i][j]*e[j]));}a[m]=p0[i]-p[i];a[(++m)++]=Q0[i]-Q[i];/*a数组存放pQ节点的不平衡量*/ }if(type[i]==2){for(j=1;j<=bus;j++){p[i]+=(e[i]*(real[i][j]*e[j]-image[i][j]*f[j])+f[i]*(real[i][j]*f[j]+image[i][j]*e[j]));}b[n]=p0[i]-p[i];b[(++n)++]=e0[i]*e0[i]+f0[i]*f0[i]-e[i]*e[i]-f[i]*f[i];}/*b数组存放pv节点的不平衡量*/}for(i=1;i<=2*pQ;i++)ub[i]=a[i];/*ub存放的是修正方程的不平衡列向量*/for(i=1;i<=2*pv;i++)ub[i+2*pQ]=b[i];for(i=1;i<=2*bus-2;i++)fprintf(fp2," ub[%d]=%f",i,ub[i]);/*不平衡量已经求出并已经打印*//*下边开始求雅克比矩阵*/fprintf(fp2,"\n------------------------------------------------------------------------------------");fprintf(fp2,"\n雅克比矩阵为:");for(i=1;i<=bus;i++)for(j=1;j<=bus;j++){a1[i][j]=0;b1[i][j]=0;}for(i=1;i<=bus;i++){for(j=1;j<=bus;j++)if(i!=j){a1[i][i]+=(real[i][j]*e[j]-image[i][j]*f[j]);b1[i][i]+=(real[i][j]*f[j]+image[i][j]*e[j]);}a1[i][i]+=(real[i][i]*e[i]-image[i][i]*f[i]);b1[i][i]+=(real[i][i]*f[i]+image[i][i]*e[i]);}for(i=1;i<=bus;i++)for(j=1;j<=bus;j++){if(type[i]==1){H[i][j]=(-1)*image[i][j]*e[i]+real[i][j]*f[i];N[i][j]=real[i][j]*e[i]+image[i][j]*f[i];J[i][j]=(-1)*N[i][j];L[i][j]=H[i][j];}else if(type[i]==2){H[i][j]=(-1)*image[i][j]*e[i]+real[i][j]*f[i];N[i][j]=real[i][j]*e[i]+image[i][j]*f[i];if(i==j){R1[i][j]=2*f[i];S[i][j]=2*e[i];}else if(i!=j){R1[i][j]=0;S[i][j]=0;}}if(i==j){H[i][j]=H[i][j]+b1[i][j];N[i][j]=N[i][j]+a1[i][j];J[i][j]=J[i][j]+a1[i][j];L[i][j]=L[i][j]-b1[i][j];}}num=1;for(i=1;i<=bus;i++){if(type[i]==1){t1[num]=i;num++;}}num=1;{if(type[j]==2){t2[num]=j;num++;}}for(i=1;i<=pQ;i++)d[i]=t1[i];for(i=pQ+1;i<=bus-1;i++)d[i]=t2[i-pQ];for(i=1;i<=2*pQ;i=i+2)for(j=1;j<=2*bus-2;j=j+2){x=d[(i+1)/2];y=d[(j+1)/2];B[i][j]=H[x][y];B[i][j+1]=N[x][y];B[i+1][j]=J[x][y];B[i+1][j+1]=L[x][y];}for(i=2*pQ+1;i<2*bus-1;i=i+2)for(j=1;j<2*bus-1;j=j+2){x=d[(i+1)/2];y=d[(j+1)/2];B[i][j]=H[x][y];B[i][j+1]=N[x][y];B[i+1][j]=R1[x][y];B[i+1][j+1]=S[x][y];}fprintf(fp2,"\n\n");for(i=1;i<2*bus-1;i++){for(j=1;j<2*bus-1;j++)fprintf(fp2," %f",B[i][j]);fprintf(fp2,"\n");}fprintf(fp2,"\n----------------------------------------------------------------------------");/*下面是用LU分解法求雅克比矩阵的逆矩阵的程序*/w=2*bus-2;/*雅克比矩阵的阶数*/for(j=1;j<=w;j++){if(i==j){E[i][j]=1;}elseE[i][j]=0;}/*单位矩阵已经生成*/for(i=1;i<=w;i++)for(j=1;j<=w;j++)zgjz[i][j]=B[i][j];for(i=1;i<=w;i++)for(j=w+1;j<=2*w;j++)zgjz[i][j]=E[i][j-w];for(i=1;i<=2*w;i++)for(j=1;j<=w;j++)Un[1][i]=zgjz[1][i];for(k=2;k<=w;k++)Ln[k][1]=zgjz[k][1]/Un[1][1];for(r=2;r<=w;r++){for(i=r;i<=2*w;i++){t=0;for(k=1;k<=r-1;k++)t+=Ln[r][k]*Un[k][i];Un[r][i]=zgjz[r][i]-t;}for(j=r+1;j<=w;j++){s=0;for(k=1;k<=r-1;k++)s+=Ln[j][k]*Un[k][r];Ln[j][r]=((zgjz[j][r]-s)/Un[r][r]);}}for(i=1;i<=w;i++)for(j=1;j<i;j++)Un[i][j]=0;for(i=1;i<=w;i++)for(j=1;j<=w;j++)Eu[i][j]=Un[i][j+w];for(j=w;j>=1;j--)Bn[w][j]=Eu[w][j]/Un[w][w];for(i=w-1;i>=1;i--)for(j=1;j<=w;j++){t=0;for(k=i+1;k<=w;k++)t+=Un[i][k]*Bn[k][j];for(k=i+1;k<=w;k++)Bn[i][j]=(Eu[i][j]-t)/Un[i][i];/*Bn[i][j]就是原雅克比矩阵的逆矩阵*/ }/*求逆矩阵的程序结束*/for(i=1;i<=w;i++)du[i]=0;for(i=1;i<=w;i++)for(j=1;j<=w;j++)du[i]+=Bn[i][j]*ub[j];fprintf(fp2,"\n修正量为:\n");for(i=1;i<=w;i++)fprintf(fp2," du[%d]=%.6f",i,du[i]);fprintf(fp2,"\n-------------------------------------------------------------------------\n");for(i=1;i<=pQ;i++){f[d[i]]+=du[(2*i-1)];e[d[i]]+=du[(2*i)];fprintf(fp2," f[%d]=%f e[%d]=%f\n",d[i],f[d[i]],d[i],e[d[i]]);}for(i=pQ+1;i<=bus-1;i++){f[d[i]]+=du[2*i-1];e[d[i]]+=du[2*i];fprintf(fp2," f[%d]=%f e[%d]=%f\n",d[i],f[d[i]],d[i],e[d[i]]);}dumax=fabs(du[1]);for(i=2;i<=w;i++){if(fabs(du[i])>dumax)dumax=fabs(du[i]);}circle++;}while(dumax>=0.00010 && circle<=50);fprintf(fp2,"\n-------------------------------------------------------------------------");fprintf(fp2,"\n ***************迭代过程结束*****************\n");fprintf(fp2,"============================================================= ===========================");fprintf(fp2,"\n平衡节点的功率为:\n");ps=0;Qs=0;for(i=1;i<=bus;i++){if(type[i]==3){for(j=1;j<=bus;j++){ps+=(real[i][j]*e[j]-image[i][j]*f[j]);Qs+=((-1)*real[i][j]*f[j]-image[i][j]*e[j]);}ps=ps*e[i];Qs=Qs*e[i];p0[i]=ps;Q0[i]=Qs;fprintf(fp2," S=%f+j%f\n",ps,Qs);}}fprintf(fp2,"pv节点的无功功率为:\n");for(i=1;i<=bus;i++){if(type[i]==2){ Q0[i]=0;for(j=1;j<=bus;j++)Q0[i]+=(f[i]*(real[i][j]*e[j]-image[i][j]*f[j])-e[i]*(real[i][j]*f[j]+image[i][j]*e[j]));fprintf(fp2," Q0[%d]=%f\n",i,Q0[i]);}}/*下边是考虑无功越限情况下的程序*/for(i=1;i<=bus;i++){if(type[i]==2 && Q0[i]>Qmax){ fprintf(fp2,"\n************************无功越限情况下的重新迭代过程*********************\n");pQ++;pv--;Q0[i]=Qmax;/*如果越限，则给pv节点的无功功率赋值为给定的Qmax*/type[i]=1;for(j=1;j<=bus;j++){p[i]=0;Q[i]=0;}goto loop;/*如果pv节点的无功功率超出了给定的Qmax，那么就返回到loop 重新进行迭代过程*/}}/*无功越限程序考虑完毕*//*下边是求线路上的功率的一段程序*/fprintf(fp2,"线路上的功率为:\n");for(i=1;i<=line;i++){real[head[i]][0]=0;real[end[i]][0]=0;image[end[i]][0]=duidi[i];image[head[i]][0]=duidi[i];}for(i=1;i<=line;i++){if(flag[i]==1)/*这是考虑单回线运行时支路的功率和损耗情况*/{pline[head[i]][end[i]]=((real[head[i]][0]+real[head[i]][end[i]])*(e[head[i]]*e[head[i]]+f[head[i]]*f [head[i]])-(e[head[i]]*e[end[i]]+f[head[i]]*f[end[i]])*real[head[i]][end[i]]-(f[head[i]]*e[end[i]]-e[ head[i]]*f[end[i]])*image[head[i]][end[i]])*(-1);Qline[head[i]][end[i]]=((-1)*(image[head[i]][0]+image[head[i]][end[i]])*(e[head[i]]*e[head[i]]+f[ head[i]]*f[head[i]])+(e[head[i]]*e[end[i]]+f[head[i]]*f[end[i]])*image[head[i]][end[i]]+(e[head[i] ]*f[end[i]]-f[head[i]]*e[end[i]])*real[head[i]][end[i]])*(-1);pline[end[i]][head[i]]=((real[end[i]][0]+real[end[i]][head[i]])*(e[end[i]]*e[end[i]]+f[end[i]]*f[end [i]])-(e[head[i]]*e[end[i]]+f[head[i]]*f[end[i]])*real[end[i]][head[i]]-(f[end[i]]*e[head[i]]-e[end[i] ]*f[head[i]])*image[end[i]][head[i]])*(-1);Qline[end[i]][head[i]]=((-1)*(image[end[i]][0]+image[end[i]][head[i]])*(e[end[i]]*e[end[i]]+f[en d[i]]*f[end[i]])+(e[head[i]]*e[end[i]]+f[head[i]]*f[end[i]])*image[end[i]][head[i]]-(f[end[i]]*e[he ad[i]]+e[end[i]]*f[head[i]])*real[end[i]][head[i]])*(-1);dpline[head[i]][end[i]]=pline[head[i]][end[i]]+pline[end[i]][head[i]];dQline[head[i]][end[i]]=Qline[head[i]][end[i]]+Qline[head[i]][end[i]];fprintf(fp2," S[%d][%d]=%f +j %f\n",head[i],end[i],pline[head[i]][end[i]],Qline[head[i]][end[i]]);fprintf(fp2," S[%d][%d]=%f +j %f\n",end[i],head[i],pline[end[i]][head[i]],Qline[end[i]][head[i]]);}else if(flag[i]==2)/*这是考虑双回线运行时支路的功率和损耗情况*/{pline[head[i]][end[i]]=pline[head[i]][end[i]]/2;Qline[head[i]][end[i]]=Qline[head[i]][end[i]]/2;pline[end[i]][head[i]]=pline[end[i]][head[i]]/2;Qline[end[i]][head[i]]=Qline[end[i]][head[i]]/2;dpline[head[i]][end[i]]=dpline[head[i]][end[i]]/2;dQline[head[i]][end[i]]=dQline[head[i]][end[i]]/2;fprintf(fp2," S[%d][%d]=%f +j %f\n",head[i],end[i],pline[head[i]][end[i]],Qline[head[i]][end[i]]);fprintf(fp2," S[%d][%d]=%f +j %f\n",end[i],head[i],pline[end[i]][head[i]],Qline[end[i]][head[i]]);}}fprintf(fp2,"\n-------------------------------------------------------------------------------------------\n");fprintf(fp2,"考虑平行之路和接地之路情况时\n单条线路损耗的功率为:\n");for(i=1;i<=line;i++)fprintf(fp2,"\n 线路%d 上的功率损耗为: %f +j %f",i,dpline[head[i]][end[i]],dQline[head[i]][end[i]]);fprintf(fp2,"\n-------------------------------------------------------------------------------------------\n");/*下边是求网络总损耗的程序*/fprintf(fp2,"网络总损耗为:\n");dpsum=0;dQsum=0;for(i=1;i<=bus;i++)dpsum+=p0[i];for(i=1;i<=bus;i++)dQsum+=Q0[i];fprintf(fp2," ds= %f +j %f\n",dpsum,dQsum);/*网络总损耗完毕*/fprintf(fp2,"-----------------------------------------------------------------------------------------------\n");/*下边是求数点效率的程序*/fprintf(fp2,"输电效率为:\n");p1=0;p2=0;for(i=1;i<=bus;i++){if(type[i]!=3){if (p0[i]>0)p1+=p0[i];elsep2+=(-1)*p0[i];}}p1=p1+ps;efficiency=p2/p1*100;/*输电效率完毕*/fprintf(fp2,"efficiency=%.4f％",efficiency);/*efficiency是输电效率*/fprintf(fp2,"\n============================================================ ===================================\n");fprintf(fp2," **********************潮流计算结束*********************");fclose(fp2);/*关闭output.txt文件*/printf("!!******计算结果存放于out.txt文件中******!!\n");}。
电力系统潮流计算C程序精编版
电力系统潮流计算注：这是一个基于N－R法的潮流计算通用程序，仅提供了子程序，需要做些处理才能成为一个可运行的计算程序！此程序非我原创，仅与大家共享！！！/******************************************************************* * 这里提供的是电力系统潮流计算机解法的五个子程序,采用的方法是 ** Newton_Raphson法.** 程序中所用的变量说明如下:** N:网络节点总数. M:网络的PQ节点数. ** L:网络的支路总数. N0:雅可比矩阵的行数. ** N1:N0+1K:打印开关.K=1,则打印;否则,不打印.** K1:子程序PLSC中判断输入电压的形式.K1=1,则为极座标形式.否则** 为直角坐标形式.** D:有功及无功功率误差的最大值. * * G(I,J):Ybus的电导元素(实部).** B(I,J):Ybus的电纳元素(虚部).** G1(I) :第I支路的串联电导. B1(I):第I支路的串联电纳. ** C1(I) :第I支路的pie型对称接地电纳. ** C(I,J):第I节点J支路不对称接地电纳. ** CO(I) :第I节点的接地电纳.** S1(I) :第I节点的起始节点号. E1(I):第I节点的终止节点号. ** P(I) :第I节点的注入有功功率. Q(I):第I节点的注入无功功率.** P0(I) :第I节点有功功率误差. Q0(I):第I节点无功功率误差. ** V0(I) :第I节点(PV节点)的电压误差(平方误差). ** V(I) :第I节点的电压误差幅值. * * E(I) :第I节点的电压的实部. F(I):第I节点的电压的虚部. ** JM(I,J):Jacoby矩阵的第I行J列元素. ** A(I,J):修正方程的增广矩阵,三角化矩阵的第I行J列元素,运算结 ** 束后A矩阵的最后一列存放修正的解. ** P1(I) :第I支路由S1(I)节点注入的有功功率. ** Q1(I) :第I支路由S1(I)节点注入的无功功率. ** P2(I) :第I支路由E1(I)节点注入的有功功率. ** Q2(I) :第I支路由E1(I)节点注入的无功功率. ** P3(I) :第I支路的有功功率损耗. * * Q3(I) :第I支路的无功功率损耗. * * ANGLE(I):第I节点电压的角度.********************************************************************/ #include <math.h>#include <stdio.h>#define f1(i) (i-1)/* 把习惯的一阶矩阵的下标转化为C语言数组下标*/#define f2(i,j,n) ((i-1)*(n)+j-1)/* 把习惯的二阶矩阵的下标转化为C语言数组下标*//***************************************************** 本子程序根据所给的支路导纳及有关信息,形成结点 ** 导纳矩阵,如打印参数K=1,则输出电导矩阵G和电纳矩B *****************************************************/void ybus(int n,int l,int m,float *g,float *b,float *g1,float *b1,float *c1,\float *c,float *co,int k,int *s1,int *e1) {extern FILE *file4;FILE *fp;int i,j,io,i0;int pos1,pos2;int st,en;if(file4==NULL){fp=stdout;}else{fp=file4; /* 输出到文件 */}/* 初始化矩阵G,B */for(i=1;i<=n;i++){for(j=1;j<=n;j++){pos2=f2(i,j,n);g[pos2]=0;b[pos2]=0;}}/* 计算支路导纳 */for(i=1;i<=l;i++){/* 计算对角元 */pos1=f1(i);st=s1[pos1];en=e1[pos1];pos2=f2(st,st,n);g[pos2]+=g1[pos1];b[pos2]+=b1[pos1]+c1[pos1];pos2=f2(en,en,n);g[pos2]+=g1[pos1];b[pos2]+=b1[pos1]+c1[pos1];/* 计算非对角元 */pos2=f2(st,en,n);g[pos2]-=g1[pos1];b[pos2]-=b1[pos1];g[f2(en,st,n)]=g[f2(st,en,n)];b[f2(en,st,n)]=b[f2(st,en,n)];}/* 计算接地支路导纳 */for(i=1;i<=n;i++){/* 对称部分*/b[f2(i,i,n)]+=co[f1(i)];/* 非对称部分*/for(j=1;j<=l;j++){b[f2(i,i,n)]+=c[f2(i,j,l)];}}if(k!=1){return; /* 如果K不为 1,则返回;否则,打印导纳矩阵 */ }fprintf(fp,"\n BUS ADMITTANCE MATRIXY(BUS):");fprintf(fp,"\n ******************* ARRAY G ********************"); for(io=1;io<=n;io+=5){i0=(io+4)>n?n:(io+4);fprintf(fp,"\n");for(j=io;j<=i0;j++){fprintf(fp,"%13d",j);}for(i=1;i<=n;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%13.6f",g[f2(i,j,n)]);}}fprintf(fp,"\n");}fprintf(fp,"\n ******************* ARRAY B ********************");for(io=1;io<=n;io+=5){i0=(io+4)>n?n:(io+4);fprintf(fp,"\n");for(j=io;j<=i0;j++){fprintf(fp,"%13d",j);}for(i=1;i<=n;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%13.6f",b[f2(i,j,n)]);}}fprintf(fp,"\n");}fprintf(fp,"\n************************************************");}/******************************************** 本子程序根据所给的功率及电压等数据 ** 求出功率及电压误差量,并返回最大有功功率 ** 以用于与给定误差比较.如打印参数K=1,则输 ** 出P0,Q0(对PQ结点),V0(对PV结点). ********************************************/void dpqc(float *p,float *q,float *p0,float *q0,float *v,float *v0,int m,\int n,float *e,float *f,int k,float *g,float *b,float *dd){extern FILE *file4;FILE *fp;int i,j,l;int pos1,pos2;float a1,a2,d1,d;if(file4==NULL){fp=stdout; /* 输出到屏幕 */}else{fp=file4; /* 输出到文件*/}l=n-1;if(k==1){fprintf(fp,"\n CHANGE OFP0,V**2,P0(I),Q0(I),V0(I) ");fprintf(fp,"\n I P0(I)Q0(I)");}for(i=1;i<=l;i++){a1=0;a2=0;pos1=f1(i);for(j=1;j<=n;j++){/* a1,a2对应课本p185式(4-86)中括号内的式子 */pos2=f2(i,j,n);a1+=g[pos2]*e[f1(j)]-b[pos2]*f[f1(j)];a2+=g[pos2]*f[f1(j)]+b[pos2]*e[f1(j)];}/* 计算式(4-86)(4-87)中的deltaPi */p0[pos1]=p[pos1]-e[pos1]*a1-f[pos1]*a2;if(i <= m){ /* 计算PQ结点中的deltaQi */q0[pos1]=q[pos1]-f[pos1]*a1+e[pos1]*a2;}else{ /* 计算PV结点中的deltaVi平方*/v0[pos1]=v[pos1]*v[pos1]-e[pos1]*e[pos1]-f[pos1]*f[pos1];}/* 输出结果 */if(k==1){if(i<m){fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]);}else if(i==m){fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]);fprintf(fp,"\n I P0(I)V0(I)");}else{fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],v0[pos1]);}}}/* 找到deltaP和deltaQ中的最小者，作为收敛指标, 存在dd中 */d=0;for(i=1;i<=l;i++){pos1=f1(i);d1=p0[pos1]>0 ? p0[pos1] : -p0[pos1];if(d<d1){d=d1;}if(i<=m){d1=q0[pos1]>0?q0[pos1]:-q0[pos1];if(d<d1){d=d1;}}}(*dd)=d;}/**************************************************** 本子程序根据节点导纳及电压求Jacoby矩阵,用于求** 电压修正量,如打印参数K=1,则输出Jacoby矩阵. ** 对应于课本P186式(4-89)(4-90) ****************************************************/void jmcc(int m,int n,int n0,float *e,float *f,float *g,float *b,float *jm,int k){extern FILE *file4;FILE *fp;int i,j,i1,io,i0,ns;int pos1,pos2;if(file4==NULL){fp=stdout;}else{fp=file4;}/* 初始化矩阵jm */for(i=1;i<=n0;i++){for(j=1;j<=n0;j++){jm[f2(i,j,n0)]=0;}}ns=n-1; /* 去掉一个平衡结点 *//* 计算式(4-89)(4-90) */for(i=1;i<=ns;i++){/* 计算式(4-90) */for(i1=1;i1<=n;i1++){/* pos1是式(4-90)中的j */pos1=f1(i1);/* pos2是式(4-90)中的ij */pos2=f2(i,i1,n);if(i<=m) /* i是PQ结点 */{/* 计算式(4-90)中的Jii等式右侧第一部分 */jm[f2(2*i-1,2*i-1,n0)]+=g[pos2]*f[pos1]+b[pos2]*e[pos1 ];/* 计算式(4-90)中的Lii等式右侧第一部分 */jm[f2(2*i-1,2*i,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1];}/* 计算式(4-90)中的Hii等式右侧第一部分 */jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1];/* 计算式(4-90)中的Nii等式右侧第一部分 */jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]-b[pos2]*e[pos1];}/* pos2是式(4-90)中的ii */pos2=f2(i,i,n);/* pos1是式(4-90)中的i */pos1=f1(i);if(i<=m) /* i是PQ结点 */{/* 计算式(4-90)中的Jii */jm[f2(2*i-1,2*i-1,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1] ;/* 计算式(4-90)中的Lii */jm[f2(2*i-1,2*i,n0)]+=g[pos2]*e[pos1]+b[pos2]*f[po s1];}/* 计算式(4-90)中的Hii */jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]-b[pos2]*f[pos1];/* 计算式(4-90)中的Jii */jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1];if(i>m) /* PV结点 */{/* 计算式(4-90)中的Rii */jm[f2(2*i-1,2*i-1,n0)]=-2*e[pos1];/* 计算式(4-90)中的Sii */jm[f2(2*i-1,2*i,n0)]=-2*f[pos1];}/* 计算式(4-89) */for(j=1;j<=ns;j++){if(j!=i){/* pos1是式(4-89)中的i */pos1=f1(i);/* pos2是式(4-89)中的ij */pos2=f2(i,j,n);/* 计算式(4-89)中的Nij */jm[f2(2*i,2*j,n0)]=b[pos2]*e[pos1]-g[pos2]*f[pos1];/* 计算式(4-89)中的Hij */jm[f2(2*i,2*j-1,n0)]=-g[pos2]*e[pos1]-b[pos2]*f[pos1];if(i<=m) /* i是PQ结点 */{/* 计算式(4-89)中的Lij (=-Hij) */jm[f2(2*i-1,2*j,n0)]=-jm[f2(2*i,2*j-1,n0)];/* 计算式(4-89)中的Jij (=Nij) */jm[f2(2*i-1,2*j-1,n0)]=jm[f2(2*i,2*j,n0)];}else /* i是PV结点 */{/* 计算式(4-89)中的Rij (=0) */jm[f2(2*i-1,2*j-1,n0)]=0;/* 计算式(4-89)中的Sij (=0) */jm[f2(2*i-1,2*j,n0)]=0;}}}}if(k!=1){return;}/* 输出Jacoby矩阵 */fprintf(fp,"\n J MATRIX(C)"); for(io=1;io<=n0;io+=5){i1=(io+4)>n0?n0:(io+4);fprintf(fp,"\n");for(j=io;j<=i1;j++){fprintf(fp,"%10d",j);}for(i=1;i<=n0;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i1;j++){fprintf(fp,"%12.6f",jm[f2(i,j,n0)]);}}}fprintf(fp,"\n");}/*********************************************** 本子程序用选列主元素的高斯消元法求解组 ** 性方程组求各结点电压修正量,如打印参数K=1,则** 输出增广矩阵变换中的上三角及电压修正量.如果** 无唯一解,则给出信息,并停止程序运行. ***********************************************/void sevc ( float a[], int n0, int k, int n1){extern FILE *file4;FILE *fp;int i,j,l,n2,n3,n4,i0,io,j1,i1;float t0,t,c;if(file4==NULL) fp=stdout;else fp=file4;for(i=1;i<=n0;i++){l=i;for(j=i;j<=n0;j++){if( fabs(a[f2(j,i,n1)]) > fabs(a[f2(l,i,n1)]) ){l=j; /* 找到这行中的最大元 */}}if(l!=i){ /* 行列交换 */for (j=i;j<=n1;j++){t=a[f2(i,j,n1)];a[f2(i,j,n1)]=a[f2(l,j,n1)];a[f2(l,j,n1)]=t;}}if (fabs(a[f2(i,i,n1)]-0)<1e-10){ /* 对角元近似于0, 无解 */printf("\nNo Solution\n");exit (1);}t0=a[f2(i,i,n1)];for(j=i;j<=n1;j++){/* 除对角元 */a[f2(i,j,n1)]/=t0;}if(i==n0){ /* 最后一行，不用消元 */continue;}/* 消元 */j1=i+1;for(i1=j1;i1<=n0;i1++){c=a[f2(i1,i,n1)];for(j=i;j<=n1;j++){a[f2(i1,j,n1)] -= a[f2(i,j,n1)] *c;}}}if(k==1){ /* 输出上三角矩阵 */fprintf(fp,"\nTrianglar Angmentex Matrix ");for(io=1;io<=n1;io+=5){i0=(io+4)>n1?n1:(io+4);fprintf(fp,"\n");fprintf(fp," ");for(i=io;i<=i0;i++){fprintf(fp,"%12d",i);}for(i=1;i<=n0;i++){fprintf(fp,"\n");fprintf(fp,"%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%15.6f", a[f2(i,j,n1)]);}}}}/* 回代求方程解 */n2=n1-2;for(i=1;i<=n2;i++){n3=n1-i;for(i1=n3;i1<=n0;i1++){n4=n0-i;a[f2(n4,n1,n1)] -= a[f2(i1,n1,n1)]*a[f2(n4,i1,n1)];}}if(k!=1){return;}/* 输出电压修正值 */fprintf(fp,"\nVoltage correction E(i), F(i) :");for(io=1;io<=n0;io+=4){i1=(io+1)/2;i0=((io+3)/2)>(n0/2)?(n0/2):((io+3)/2);fprintf(fp,"\n");for(j=i1;j<=i0;j++){fprintf(fp,"%16d%16d",j,j);}i1 = 2*i0;fprintf(fp,"\n");for(i=io;i<=i1;i++){fprintf(fp,"%15.6f", a[f2(i,n1,n1)]);}}}#define Pi 3.1415927/180void plsc(int n,int l,int m,float g[],float b[],float e[],float f[],\int e1[],int s1[],float g1[],floatb1[],float c1[],float c[],\float co[],float p1[],float q1[],floatp2[],float q2[],float p3[],\float q3[],float p[],float q[],floatv[],float angle[],int k1){extern FILE *file4;FILE *fp;float t1,t2,st,en,cm,x,y,z,x1,x2,y1,y2;int i,i1,j,m1,ns,pos1,pos2,km;ns=n-1;if(file4==NULL){fp=stdout;}else{fp=file4;}fprintf(fp,"\nTHE RESULT ARE:");if(k1==1){for(i=0;i<n;i++){angle[i]*=Pi;e[i]=v[i]*cos(angle[i]);f[i]=v[i]*sin(angle[i]);}}t1=0.0;t2=0.0;for(i=1;i<=n;i++){pos1=f1(i);pos2=f2(n,i,n);t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];}pos1=f1(n);p[pos1]=t1*e[pos1];q[pos1]=-t2*e[pos1];m1=m+1;for(i1=m1;i1<=ns;i1++){t1=0;t2=0;for(i=1;i<=n;i++){pos1=f1(i);pos2=f2(i1,i,n);t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];}pos1=f1(i1);q[pos1]=f[pos1]*t1-e[pos1]*t2;}for(i=0;i<n; i++){cm=co[i];if(cm!=0){q[i]-=(e[i]*e[i]+f[i]*f[i])*cm;}}fprintf(fp,"\nBUS DATA");fprintf(fp,"\nBUS VOLTAGE ANGLE(DEG S.) BUS P BUS Q");for(i=0;i<n;i++){v[i]=sqrt(e[i]*e[i]+f[i]*f[i]);x=e[i];y=f[i];z=y/x;angle[i]=atan(z);angle[i]/=Pi;fprintf(fp,"\n%3d%13.5e%15.5f%15.5e%15.5e",i+1,v[i],angle[ i],p[i],q[i]);}fprintf(fp,"\n LINE FLOW ");for(i=1;i<=l;i++){pos1=f1(i);st=s1[pos1];en=e1[pos1];x1=e[f1(st)]*e[f1(st)]+f[f1(st)]*f[f1(st)];x2=e[f1(en)]*e[f1(en)]+f[f1(en)]*f[f1(en)];y1=e[f1(st)]*e[f1(en)]+f[f1(st)]*f[f1(en)];y2=f[f1(st)]*e[f1(en)]-e[f1(st)]*f[f1(en)];p1[pos1]=(x1-y1)*g1[pos1]-y2*b1[pos1];q1[pos1]=-x1*(c1[pos1]+b1[pos1])+y1*b1[pos1]-y2*g1[pos1];p2[pos1]=(x2-y1)*g1[pos1]+y2*b1[pos1];q2[pos1]=-x2*(c1[pos1]+b1[pos1])+y1*b1[pos1]+y2*g1[pos1];for(j=1;j<=n;j++){cm=c[f2(j,i,l)];if(cm!=0.0){km=1;if(en==j){km=2;}if(km==1){q1[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;}else{q2[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;}}}p3[pos1]=p1[pos1]+p2[pos1] ;q3[pos1]=q1[pos1]+q2[pos1] ;fprintf(fp,"\n%2d%8d%11d%13.6e%13.6e%13.6e%13.6e%17d%11d%1 3.6e%13.6e",\i,s1[pos1],e1[pos1],p1[pos1],q1[pos1],p3[pos1],q3[pos1 ],\e1[pos1],s1[pos1],p2[pos1],q2[pos1]);}}。
电力系统潮流计算程序
电力系统潮流计算c语言程序，两行，大家可以看看，仔细研究，然后在这个基础上修改。
谢谢#include "stdafx.h"#include <iostream>#include <fstream>#include <process.h>#include"Complex.h"#include"wanjing.h"#include"gauss.h"using namespace std;int _tmain(int argc, _TCHAR* argv[]){int i; //i作为整个程序的循环变量int N=Bus::ScanfBusNo(); //输入节点个数int L=Line::ScanflineNo(); //输入支路个数if((L&&N)==0){return 0;} //如果找不到两个文件中的任意一个，退出Line *line=new Line[L]; //动态分配支路结构体Line::ScanfLineData(line); //输入支路参数Line::PrintfLineData(line,L); //输出支路参数Bus *bus=new Bus[N]; //动态分配结点结构体for(int i=0;i<N;i++){bus[i].Sdelta.real=0;bus[i].Sdelta.image=0;}Bus::ScanfBusData(bus); //输入节点参数Bus::PrintfBusData(bus,N); //输出结点参数Complex **X; X=new Complex *[N]; for(i=0;i<N;i++) {X[i]=new Complex[N];}//动态分配二维数组作节点导纳阵Bus::JisuanNodeDnz(X,line,bus,L,N); //计算节点导纳矩阵Bus::PrintfNodeDnz(X,N); //输出节点导纳矩阵int NN=(N-1)*2;double **JacAug; JacAug=new double *[NN]; for(i=0;i<NN;i++) {JacAug[i]=new double[NN+1];}//为雅可比增广动态分配一个二维数组double *x;x=new double[NN];int count=1;LOOP:Bus::JisuanNodeI(X,bus,N); //计算节点注入电流Bus::JisuanNodeScal(X,bus,N); //计算节点功率Bus::JisuanNodeScal(X,bus,N);//计算节点功率Bus::JisuanNodeSdelta(bus,N); //计算节点功率差值Bus::PrintfNodeScal(X,bus,N);//输出节点功率差值int icon=wehcon1(bus,N);//whether converbence看迭代是否结束if(icon==1){cout<<"icon="<<icon<<"进行第"<<count<<"次迭代"<<endl;Bus::JisuanJacAug(JacAug,X,bus,N); //计算雅可比增广矩阵// Bus::PrintfJacAug(JacAug,N);gauss::gauss_slove(JacAug,x,NN);//解方程组求出电压差值Bus::ReviseNodeV(bus,x,N);//修正节点电压// Bus::PrintfNodeV(bus,N);count++;goto LOOP;}else{for(i=0;i<L;i++)//计算输电线路上的传输功率{int statemp,endtemp;Complex aa,bb,cc,dd,B;B.real=0;B.image=-line[i].B;statemp=line[i].start;endtemp=line[i].end;aa=Complex::productComplex (Complex::getconj(bus[statemp-1].V), B);bb=Complex::subComplex (Complex::getconj(bus[statemp-1].V), Complex::getconj(bus[endtemp-1].V));cc=Complex::productComplex (bb , Complex::getconj(line[i].Y));dd=Complex::CaddC(aa,cc);line[i].stoe=Complex::productComplex(b us[statemp-1].V,dd);aa=Complex::productComplex (Complex::getconj(bus[endtemp-1].V), B);bb=Complex::subComplex (Complex::getconj(bus[endtemp-1].V), Complex::getconj(bus[statemp-1].V));cc=Complex::productComplex (bb , Complex::getconj(line[i].Y));dd=Complex::CaddC(aa,cc);line[i].etos=Complex::productComplex(b us[endtemp-1].V,dd);}cout<<"icon="<<icon<<",迭代结束。
电力系统通用潮流计算C语言程序
#include <iostream>#include <fstream>#include<iomanip>#include<math.h>using namespace std;//节点号类型负荷有功负荷无功母线数据（类型1＝PV节点，2＝PQ节点，3＝平衡节点）struct BUS{int busno;int type;float Pd;float Qd;};//发电机数据节点号有功发电电压幅值struct Generator{int busno;float Pg;float Vg;};//支路信息节点I 节点J R X B/2 kstruct Line{int busi;int busj;float R;float X;float B;float k;};//deltaP deltaQ deltaV^2//void fun1(double YG[][50],double YB[][50],double e[],double f[],int type[],int N,double W[],double P[],double Q[],double V[]){double dP=0,dQ=0,dV=0;int i,j;for(i=0;i<N-1;i++){double A=0,B=0;for(j=0;j<N;j++){A+=YG[i][j]*e[j]-YB[i][j]*f[j];B+=YG[i][j]*f[j]+YB[i][j]*e[j];}dV=V[i]*V[i]-e[i]*e[i]-f[i]*f[i];dP=P[i]-e[i]*A-f[i]*B;W[2*i]=dP;dQ=Q[i]-f[i]*A+e[i]*B;if(type[i]==1)W[2*i+1]=dQ;else W[2*i+1]=dV;}}//Jacobi矩阵//void Jacobi(double YG[][50],double YB[][50],double e[50],double f[50],int type[50],int N ,double Ja[100][101]){ int i,j;for(i=0;i<N;i++){for(j=0;j<N;j++){if(i!=j){if(type[i]==1){Ja[2*i][2*j]=-(YG[i][j]*e[i]+YB[i][j]*f[i]);Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i];Ja[2*i+1][2*j]=Ja[2*i][2*j+1];Ja[2*i+1][2*j+1]=-Ja[2*i][2*j];}else {Ja[2*i][2*j]=-YG[i][j]*e[i]+YB[i][j]*f[i];Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i];Ja[2*i+1][2*j+1]=Ja[2*i+1][2*j]=0;}}else {double a[50]={0},b[50]={0};for(int k=0;k<N;k++){a[i]+=(YG[i][k]*e[k]-YB[i][k]*f[k]);b[i]+=(YG[i][k]*f[k]+YB[i][k]*e[k]);Ja[2*i][2*j]=-a[i]-YG[i][i]*e[i]-YB[i][i]*f[i];Ja[2*i][2*j+1]=-b[i]+YB[i][i]*e[i]-YG[i][i]*f[i];if(type[i]==1){Ja[2*i+1][2*j]=b[i]+YB[i][i]*e[i]-YG[i][i]*f[i];Ja[2*i+1][2*j+1]=-a[i]+YG[i][i]*e[i]+YB[i][i]*f[i];}else {Ja[2*i+1][2*j]=-2*e[i];Ja[2*i+1][2*j+1]=-2*f[i];}}}}}}//高斯消元法解方程组函数//void gauss(double a[][101],int n){int i,j,k;double c;for(k=0;k<n-1;k++) {c=a[k][k];for(j=k;j<=n;j++) a[k][j]/=c;for(i=k+1;i<n;i++) {c=a[i][k];for(j=k;j<=n;j++) a[i][j]-=c*a[k][j];}}a[n-1][n]/=a[n-1][n-1];for(k=n-2;k>=0;k--)for(j=k+1;j<n;j++) a[k][n]-=a[k][j]*a[j][n];}void main(){ifstream fin;int N=0,GS=0,LD=0,ZLs=0; //节点数发电机数负荷数支路数// BUS *B;Generator *G;Line *L;//从文本中读入原始数据到数组中//fin.open("C:\\data.txt");if(!fin){cout<<"输入数据文件不存在！"<<endl;getchar();}int m1[50]={0},m2[50]={0};float m3[50],m4[50],m5[50],m6[50];int i,j,l;for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m2[i]>>m3[i]>>m4[i];N++;}B =new BUS[N];for (i=0;i<N;i++){B[i].busno=m1[i];B[i].type=m2[i];B[i].Pd=m3[i];B[i].Qd=m4[i];}for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m4[i]>>m3[i];GS++;}G =new Generator[GS];for (i=0;i<GS;i++){G[i].busno=m1[i];G[i].Pg=m4[i];G[i].Vg=m3[i];}for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m2[i]>>m3[i]>>m4[i]>>m5[i]>>m6[i];ZLs++;}L =new Line[ZLs];for (i=0;i<ZLs;i++){L[i].busi=m1[i];L[i].busj=m2[i];L[i].R=m3[i];L[i].X=m4[i];L[i].B=m5[i];L[i].k=m6[i];}LD=N-GS;fin.close();//节点导纳矩阵形成//double YB[50][50],YG[50][50],BB[50][50],K[50][50];for(i=0;i<N;i++){for(j=0;j<N;j++){YB[i][j]=0;YG[i][j]=0;BB[i][j]=0;K[i][j]=1;}}for (l=0;l<ZLs;l++){i=L[l].busi-1;j=L[l].busj-1;K[i][j]=L[l].k;BB[i][j]=BB[j][i]=L[l].B;YG[i][j]=YG[j][i]=L[l].R/(L[l].R*L[l].R+L[l].X*L[l].X);YB[i][j]=YB[j][i]=-L[l].X/(L[l].R*L[l].R+L[l].X*L[l].X);}for(i=0;i<N;i++){for(j=i;j<N;j++) {K[i][j]=K[j][i];K[j][i]=1;}for(j=0;j<N;j++){if(i!=j){YG[i][i]=YG[i][i]+(YG[i][j]*K[i][j]*K[i][j]);YB[i][i]=YB[i][i]+(YB[i][j]*K[i][j]*K[i][j]+BB[i][j]);}}}//修正后//for (l=0;l<ZLs;l++){i=L[l].busi-1;j=L[l].busj-1;K[i][j]=L[l].k;YG[i][j]=-YG[i][j]*K[i][j];YG[j][i]=YG[i][j];YB[i][j]=-YB[i][j]*K[i][j];YB[j][i]=YB[i][j];}int type[50]={0};for(i=0;i<N;i++){type[i]=B[i].type;}//PQV的获得//double P[50],Q[50],V[50];for(i=0;i<N;i++){P[i]=0;Q[i]=0;V[i]=0;P[i]=-B[i].Pd;Q[i]=-B[i].Qd;}for (i=0;i<GS;i++){P[G[i].busno-1]=G[i].Pg;V[G[i].busno-1]=G[i].Vg;}// 求A=e+f//double e[50]={0},f[50]={0};double C[100]={0},D[100]={0};for(i=0;i<N;i++){if(V[i]==0){C[2*i]=1;}else C[2*i]=V[i];}double W[100]={0},Ja[100][101]={0}; //调用Jacobi函数和高斯函数//for(int t=1;t<10;t++){for(i=0;i<2*N-2;i++){e[i]=C[2*i];f[i]=C[2*i+1];}fun1(YG,YB,e,f,type,N,W,P,Q,V);double it=fabs(W[0]);for(i=1;i<2*N-2;i++){if (it<fabs(W[i])) {it=fabs(W[i]);j=i;}}//中间迭代过程//cout<<setw(10)<<"迭代次数"<<setw(20)<<"最大的功率误差"<<setw(8)<<"节点号"<<endl;cout<<setw(10)<<t<<setw(20)<<it<<setw(8)<<j/2+1<<endl;if (it<0.00001) break;Jacobi(YG,YB,e,f,type,N,Ja);for(i=0;i<2*N-2;i++){Ja[i][2*N-2]=W[i];}//高斯消元法解方程//gauss(Ja,2*N-2);for(i=0;i<2*N-2;i++){D[i]=-Ja[i][2*(N-1)];C[i]+=D[i];}}//平衡节点//for(i=0;i<N;i++){double a=0,b=0;for(int j=0;j<N;j++){a+=(YG[i][j]*e[j]-YB[i][j]*f[j]);b+=(YB[i][j]*e[j]+YG[i][j]*f[j]);}P[i]=e[i]*a+f[i]*b;Q[i]=f[i]*a-e[i]*b;}//支路//double PZL[100][101]={0},QZL[100][101]={0},pr[100][101]={0},qx[100][101]={0};double x1=0,x2=0,y1=0,y2=0,I2=0;for(int k=0;k<ZLs;k++){i=L[k].busi-1;j=L[k].busj-1;x1=e[i]/L[k].k-e[j];y1=f[i]/L[k].k-f[j];x2=-e[i]*YG[i][j]-f[i]*YB[i][j];y2=-f[i]*YG[i][j]+e[i]*YB[i][j];QZL[i][j]=(x1*y2-x2*y1);PZL[i][j]=(x1*x2+y1*y2);I2=(PZL[i][j]*PZL[i][j]+QZL[i][j]*QZL[i][j])/(e[i]*e[i]+f[i]*f[i]);pr[i][j]=I2*L[k].R;qx[i][j]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[i][j]+=(e[i]*e[i]+f[i]*f[i])*(-L[k].B);x1=e[j]*L[k].k-e[i];y1=f[j]*L[k].k-f[i];x2=-e[j]*YG[j][i]-f[j]*YB[j][i];y2=-f[j]*YG[j][i]+e[j]*YB[j][i];QZL[j][i]=(x1*y2-x2*y1);PZL[j][i]=(x1*x2+y1*y2);I2=(PZL[j][i]*PZL[j][i]+QZL[j][i]*QZL[j][i])/(e[j]*e[j]+f[j]*f[j]);pr[j][i]=I2*L[k].R;qx[j][i]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[j][i]+=(e[j]*e[j]+f[j]*f[j])*(-L[k].B);}//全网数据//int high=1,low=1;double PG=0,PL=0,Prr=0,Vh=sqrt(e[0]*e[0]+f[0]*f[0]),Vl=sqrt(e[0]*e[0]+f[0]*f[0]);for(k=0;k<N;k++){Prr+=P[k];if(B[k].type==1) PL+=B[k].Pd;else PG+=P[k];if(sqrt(e[k]*e[k]+f[k]*f[k])>Vh){Vh=sqrt(e[k]*e[k]+f[k]*f[k]);high=k+1;}if(sqrt(e[k]*e[k]+f[k]*f[k])<Vl){Vl=sqrt(e[k]*e[k]+f[k]*f[k]);low=k+1;}}//输出数据到文件databak.txt//ofstream fout;fout.open("C:\\databak.txt");fout<<"节点"<<endl;fout<<setw(8)<<"节点号"<<setw(16)<<"V"<<setw(16)<<"弧度"<<setw(16)<<"发电P"<<setw(16)<<"发电Q"<<setw(16)<<"负荷P"<<setw(16)<<"负荷Q"<<endl;for(i=0;i<LD;i++){fout<<setw(8)<<i+1<<setw(16)<<sqrt(e[i]*e[i]+f[i]*f[i])<<setw(16)<<atan2(f[i],e[i])*180/3.14159<<setw(16)<<0 <<setw(16)<<0<<setw(16)<<B[i].Pd<<setw(16)<<B[i].Qd<<endl;}for(j=0;j<GS;j++){i=G[j].busno-1;fout<<setw(8)<<i+1<<setw(16)<<V[i]<<setw(16)<<atan2(f[i],e[i])*180/3.14159<<setw(16)<<P[i]<<setw(16)<<Q[i]<< setw(16)<<0<<setw(16)<<0<<endl;}fout<<"支路"<<endl;fout<<setw(4)<<"i"<<setw(4)<<"j"<<setw(10)<<"i_j有功"<<setw(10)<<"i_j无功"<<setw(10)<<"j_i有功"<<setw(12)<<"j_i无功"<<setw(12)<<"有功损耗"<<setw(12)<<"无功损耗"<<endl;for (k=0;k<ZLs;k++){i=L[k].busi-1;j=L[k].busj-1;fout<<setw(4)<<L[k].busi<<setw(4)<< L[k].busj <<setw(10)<<PZL[i][j] <<setw(10)<<QZL[i][j]<<setw(10)<< PZL[j][i] <<setw(12)<<QZL[j][i]<<setw(12)<<pr[i][j]<<setw(12)<<qx[i][j]<<endl;}fout<<"全网数据"<<endl;fout<<setw(14)<<"发电有功"<<setw(14)<<"负荷有功"<<setw(14)<<"有功损耗"<<setw(14)<<"最高电压"<<setw(14)<<"节点号"<<setw(14)<<"最低电压"<<setw(14)<<"节点号"<<endl;fout<<setw(14)<<PG<<setw(14)<<PL<<setw(14)<<Prr<<setw(14)<<Vh<<setw(14)<<high<<setw(14)<<Vl<<setw(14) <<low<<endl;fout.close();}。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
电力系统通用潮流计算C语言程序<fstream>#include<iomanip>#include<math、h>using namespace std; //节点号类型负荷有功负荷无功母线数据（类型1＝PV节点，2＝PQ节点，3＝平衡节点）struct BUS {int busno;int type;float Pd;float Qd;};//发电机数据节点号有功发电电压幅值struct Generator{int busno;float Pg;float Vg;};//支路信息节点I 节点J R X B/2 kstruct Line{int busi;int busj;float R;float X;float B;float k;};//deltaP deltaQ deltaV^2//void fun1(double YG[][50],doubleYB[][50],double e[],double f[],int type[],int N,doubleW[],double P[],double Q[],double V[]){doubledP=0,dQ=0,dV=0;int i,j;for(i=0;i<N-1;i++){doubleA=0,B=0;for(j=0;j<N;j++){A+=YG[i][j]*e[j]-YB[i][j]*f[j]; B+=YG[i][j]*f[j]+YB[i][j]*e[j]; } dV=V[i]*V[i]-e[i]*e[i]-f[i]*f[i]; dP=P[i]-e[i]*A-f[i]*B; W[2*i]=dP; dQ=Q[i]-f[i]*A+e[i]*B; if(type[i]==1)W[2*i+1]=dQ;else W[2*i+1]=dV;}}//Jacobi矩阵//void Jacobi(double YG[][50],double YB[][50],doublee[50],double f[50],int type[50],int N ,doubleJa[100][101]){ int i,j;for(i=0;i<N;i++){ for(j=0;j<N;j++){ if(i!=j){ if(type[i]==1){ Ja[2*i][2*j]=-(YG[i][j]*e[i]+YB[i][j]*f[i]);Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i];Ja[2*i+1][2*j]=Ja[2*i][2*j+1]; Ja[2*i+1][2*j+1]=-Ja[2*i][2*j]; } else { Ja[2*i][2*j]=-YG[i][j]*e[i]+YB[i][j]*f[i];Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i];Ja[2*i+1][2*j+1]=Ja[2*i+1][2*j]=0; } } else { doublea[50]={0},b[50]={0}; for(intk=0;k<N;k++){ a[i]+=(YG[i][k]*e[k]-YB[i][k]*f[k]);b[i]+=(YG[i][k]*f[k]+YB[i][k]*e[k]); Ja[2*i][2*j]=-a[i]-YG[i][i]*e[i]-YB[i][i]*f[i]; Ja[2*i][2*j+1]=-b[i]+YB[i][i]*e[i]-YG[i][i]*f[i];if(type[i]==1){ Ja[2*i+1][2*j]=b[i]+YB[i][i]*e[i]-YG[i][i]*f[i]; Ja[2*i+1][2*j+1]=-a[i]+YG[i][i]*e[i]+YB[i][i]*f[i]; } else{ Ja[2*i+1][2*j]=-2*e[i]; Ja[2*i+1][2*j+1]=-2*f[i]; } } } } }}//高斯消元法解方程组函数//voidgauss(double a[][101],int n){ int i,j,k;double c; for(k=0;k<n-1;k++){ c=a[k][k]; for(j=k;j<=n;j++)a[k][j]/=c; for(i=k+1;i<n;i++){ c=a[i][k]; for(j=k;j<=n;j++)a[i][j]-=c*a[k][j]; } } a[n-1][n]/=a[n-1][n-1];for(k=n-2;k>=0;k--)for(j=k+1;j<n;j++)a[k][n]-=a[k][j]*a[j][n]; } void main(){ifstream fin;int N=0,GS=0,LD=0,ZLs=0; //节点数发电机数负荷数支路数//BUS *B; Generator *G;Line *L;//从文本中读入原始数据到数组中//fin、open("C:\\data、txt");if(!fin){ cout<<"输入数据文件不存在！"<<endl;getchar();} int m1[50]={0},m2[50]={0};floatm3[50],m4[50],m5[50],m6[50];inti,j,l;for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m2[ i]>>m3[i]>>m4[i];N++;}B =new BUS[N];for(i=0;i<N;i++){B[i]、busno=m1[i];B[i]、type=m2[i];B[i]、Pd=m3[i];B[i]、Qd=m4[i];}for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m4[i]>>m3 [i];GS++;}G =new Generator[GS];for (i=0;i<GS;i++){G[i]、busno=m1[i];G[i]、Pg=m4[i];G[i]、Vg=m3[i];}for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin> >m2[i]>>m3[i]>>m4[i]>>m5[i]>>m6[i];ZLs++;}L =newLine[ZLs];for (i=0;i<ZLs;i++){L[i]、busi=m1[i];L[i]、busj=m2[i];L[i]、R=m3[i];L[i]、X=m4[i];L[i]、B=m5[i];L[i]、k=m6[i]; }LD=N-GS;fin、close();//节点导纳矩阵形成//doubleYB[50][50],YG[50][50],BB[50][50],K[50][50];for(i=0;i<N;i+ +){for(j=0;j<N;j++){YB[i][j]=0;YG[i][j]=0;BB[i][j]=0;K[i] [j]=1;}}for (l=0;l<ZLs;l++){i=L[l]、busi-1;j=L[l]、busj-1;K[i][j]=L[l]、k;BB[i][j]=BB[j][i]=L[l]、B;YG[i][j]=YG[j][i]=L[l]、R/(L[l]、R*L[l]、R+L[l]、X*L[l]、X);YB[i][j]=YB[j][i]=-L[l]、X/(L[l]、R*L[l]、R+L[l]、X*L[l]、X);}for(i=0;i<N;i++){for(j=i;j<N;j++) {K[i][j]=K[j][i];K[j][i]=1;}for(j=0;j<N;j++){if(i!=j) {YG[i][i]=YG[i][i]+(YG[i][j]*K[i][j]*K[i][j]);YB[i][i]=YB [i][i]+(YB[i][j]*K[i][j]*K[i][j]+BB[i][j]);}}}//修正后//for (l=0;l<ZLs;l++){i=L[l]、busi-1;j=L[l]、busj-1;K[i][j]=L[l]、k;YG[i][j]=-YG[i][j]*K[i][j];YG[j][i]=YG[i][j];YB[i][j]=-YB[i][j]*K[i][j];YB[j][i]=YB[i][j];}inttype[50]={0};for(i=0;i<N;i++){type[i]=B[i]、type;}//PQV的获得//doubleP[50],Q[50],V[50];for(i=0;i<N;i++){P[i]=0;Q[i]=0;V[i]=0;P [i]=-B[i]、Pd;Q[i]=-B[i]、Qd;}for (i=0;i<GS;i++){P[G[i]、busno-1]=G[i]、Pg;V[G[i]、busno-1]=G[i]、Vg;}// 求A=e+f//double e[50]={0},f[50]={0};doubleC[100]={0},D[100]={0};for(i=0;i<N;i++){if(V[i]==0){C[2*i]=1;}else C[2*i]=V[i];}doubleW[100]={0},Ja[100][101]={0};//调用Jacobi函数和高斯函数//for(int t=1;t<10;t++){for(i=0;i<2*N-2;i++){e[i]=C[2*i];f[i]=C[2*i+1];}fun1(YG,YB,e,f,type,N,W ,P,Q,V);double it=fabs(W[0]);for(i=1;i<2*N-2;i++){if(it<fabs(W[i])){it=fabs(W[i]);j=i;}}//中间迭代过程//cout<<setw(10)<<"迭代次数"<<setw(20)<<"最大的功率误差"<<setw(8)<<"节点号"<<endl;cout<<setw(10)<<t<<setw(20)<<it<<setw(8)<<j/2+1<< endl;if (it<0、00001)break; Jacobi(YG,YB,e,f,type,N,Ja);for(i=0;i<2*N-2;i++){Ja[i][2*N-2]=W[i];} //高斯消元法解方程//gauss(Ja,2*N-2); for(i=0;i<2*N-2;i++){D[i]=-Ja[i][2*(N-1)];C[i]+=D[i];}}//平衡节点//for(i=0;i<N;i++){doublea=0,b=0;for(int j=0;j<N;j++){a+=(YG[i][j]*e[j]-YB[i][j]*f[j]);b+=(YB[i][j]*e[j]+YG[i][j]*f[j]);}P[i]=e[i ]*a+f[i]*b;Q[i]=f[i]*a-e[i]*b; }//支路//doublePZL[100][101]={0},QZL[100][101]={0},pr[100][101]={0},qx[1 00][101]={0};double x1=0,x2=0,y1=0,y2=0,I2=0;for(intk=0;k<ZLs;k++){ i=L[k]、busi-1;j=L[k]、busj-1;x1=e[i]/L[k]、k-e[j];y1=f[i]/L[k]、k-f[j];x2=-e[i]*YG[i][j]-f[i]*YB[i][j];y2=-f[i]*YG[i][j]+e[i]*YB[i][j];QZL[i][j]=(x1*y2-x2*y1);PZL[i][j]=(x1*x2+y1*y2);I2=(PZL[i][j]*PZL[i][j]+QZ L[i][j]*QZL[i][j])/(e[i]*e[i]+f[i]*f[i]);pr[i][j]=I2*L[k]、R;qx[i][j]=I2*L[k]、X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k]、B;QZL[i][j]+=(e[i]*e[i]+f[i]*f[i])*(-L[k]、B);x1=e[j]*L[k]、k-e[i];y1=f[j]*L[k]、k-f[i];x2=-e[j]*YG[j][i]-f[j]*YB[j][i];y2=-f[j]*YG[j][i]+e[j]*YB[j][i];QZL[j][i]=(x1*y2-x2*y1);PZL[j][i]=(x1*x2+y1*y2);I2=(PZL[j][i]*PZL[j][i]+QZ L[j][i]*QZL[j][i])/(e[j]*e[j]+f[j]*f[j]);pr[j][i]=I2*L[k]、R;qx[j][i]=I2*L[k]、X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k]、B;QZL[j][i]+=(e[j]*e[j]+f[j]*f[j])*(-L[k]、B);} //全网数据//int high=1,low=1; doublePG=0,PL=0,Prr=0,Vh=sqrt(e[0]*e[0]+f[0]*f[0]),Vl=sqrt(e[0] *e[0]+f[0]*f[0]);for(k=0;k<N;k++){Prr+=P[k];if(B[k]、type==1)PL+=B[k]、Pd;elsePG+=P[k];if(sqrt(e[k]*e[k]+f[k]*f[k])>Vh){Vh=sqrt(e[k]*e[ k]+f[k]*f[k]);high=k+1;}if(sqrt(e[k]*e[k]+f[k]*f[k])<Vl){Vl=sqrt(e[k]*e[k]+f[k]*f[k]);low=k+1;}} //输出数据到文件databak、txt//ofstream fout;fout、open("C:\\databak、txt");fout<<"节点"<<endl; fout<<setw(8)<<"节点号"<<setw(16)<<"V"<<setw(16)<<"弧度"<<setw(16)<<"发电P"<<setw(16)<<"发电Q"<<setw(16)<<"负荷P"<<setw(16)<<"负荷Q"<<endl;for(i=0;i<LD;i++){ fout<<setw(8)<<i+1<<setw(16)< <sqrt(e[i]*e[i]+f[i]*f[i])<<setw(16)<<atan2(f[i],e[i])*18 0/3、14159<<setw(16)<<0<<setw(16)<<0<<setw(16)<<B[i]、Pd<<setw(16)<<B[i]、Qd<<endl;}for(j=0;j<GS;j++){ i=G[j]、busno-1;fout<<setw(8)<<i+1<<setw(16)<<V[i]<<setw(16)<<atan2(f[i], e[i])*180/3、14159<<setw(16)<<P[i]<<setw(16)<<Q[i]<<setw(16)<<0<<setw( 16)<<0<<endl;} fout<<"支路 "<<endl;fout<<setw(4)<<"i"<<setw(4)<<"j"<<setw(10)<<"i_j有功"<<setw(10)<<"i_j无功"<<setw(10)<<"j_i有功"<<setw(12)<<"j_i无功"<<setw(12)<<"有功损耗"<<setw(12)<<"无功损耗"<<endl; for (k=0;k<ZLs;k++){ i=L[k]、busi-1;j=L[k]、busj-1; fout<<setw(4)<<L[k]、busi<<setw(4)<<L[k]、busj <<setw(10)<<PZL[i][j]<<setw(10)<<QZL[i][j]<<setw(10)<< PZL[j][i]<<setw(12)<<QZL[j][i]<<setw(12)<<pr[i][j]<<setw(12)<<qx[i ][j]<<endl; } fout<<"全网数据"<<endl; fout<<setw(14)<<"发电有功"<<setw(14)<<"负荷有功"<<setw(14)<<"有功损耗"<<setw(14)<<"最高电压"<<setw(14)<<"节点号"<<setw(14)<<"最低电压"<<setw(14)<<"节点号"<<endl;fout<<setw(14)<<PG<<setw(14)<<PL<<setw(14)<<Prr<<setw(14) <<Vh<<setw(14)<<high<<setw(14)<<Vl<<setw(14)<<low<<endl; fout、close();}。