基于MATLAB的电力系统潮流计算%简单潮流计算的小程序,相关的原始数据数据数据输入格式如下:%B1是支路参数矩阵,第一列和第二列是节点编号。
节点编号由小到大编写%对于含有变压器的支路,第一列为低压侧节点编号,第二列为高压侧节点%编号,将变压器的串联阻抗置于低压侧处理。
%第三列为支路的串列阻抗参数。
%第四列为支路的对地导纳参数。
%第五烈为含变压器支路的变压器的变比%第六列为变压器是否是否含有变压器的参数,其中“1”为含有变压器,%“0”为不含有变压器。
%B2为节点参数矩阵,其中第一列为节点注入发电功率参数;第二列为节点%负荷功率参数;第三列为节点电压参数;第六列为节点类型参数,其中%“1”为平衡节点,“2”为PQ节点,“3”为PV节点参数。
%X为节点号和对地参数矩阵。
其中第一列为节点编号,第二列为节点对地%参数。
n=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');X=input('节点号和对地参数:X=');Y=zeros(n);Times=1; %置迭代次数为初始值%创建节点导纳矩阵for i=1:n1if 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));endendYOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %将OrgS、DetaS初始化%创建OrgS,用于存储初始功率参数h=0;j=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*r eal(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))* imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real( B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag (B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:n %对PV节点的处理,注意这时不可再将h初始化为0if i~=isb&B2(i,6)==3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*r eal(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))* imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real( B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag (B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendOrgS%创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:nif B2(i,6)==3t=t+1;PVU(t,1)=B2(i,3);endendPVU%创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);endendt=0;for i=1:n %对PV节点的处理,注意这时不可再将h初始化为0if i~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^ 2-imag(B2(i,3))^2;endendDetaS%创建I,用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3)); endendI%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理if B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1)); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1)); else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1) %将用于内循环的指针置于初始值,以确保雅可比矩阵换行k=0;endendendendendk=0;for i=1:n %对PV节点的处理if B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1) %将用于内循环的指针置于初始值,以确保雅可比矩阵换行k=0;endendendendendJacbi%求解修正方程,获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:n %对PQ节点处理if B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:n %对PV节点的处理if B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendB2%开始循环********************************************************************* *while abs(max(DetaU))>prOrgS=zeros(2*n-2,1); %!!!初始功率参数在迭代过程中是不累加的,所以在这里必须将其初始化为零矩阵h=0;j=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*r eal(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))* imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real( B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag (B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:nif i~=isb&B2(i,6)==3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*r eal(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))* imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real( B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag (B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendOrgS%创建DetaSh=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);endendt=0;for i=1:nif i~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^ 2-imag(B2(i,3))^2;endendDetaS%创建Ii=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));endendI%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1)); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1)); elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+re al(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag( Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+re al(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag( Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendendendJacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:nif B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:nif B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendB2Times=Times+1; %迭代次数加1endTimes一个原始数据的例子节点数 5支路数 5平衡节点编号 5精度pr 0.000001B1(支路参数矩阵)[1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;42 0.015i 0 1.05 1;53 0.03i 0 1.05 1]B2(节点参数矩阵)[0 -1.6-0.8i 1 0 0 2;0 -2-1i 1 0 0 2;0 -3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1]X(节点号和对地参数)[1 0;2 0;3 0;4 0;5 0]电力系统潮流计算——9结点算例-PQ法原始数据录入data.txt文档:标号,起始结点,终止结点,支路电阻参数,支路电抗参数,支路对地导纳参数1,2,5,0.0,0.063,0.0,2,5,9,0.019,0.072,0.075,3,6,9,0.012,0.101,0.105,4,3,6,0.0,0.059,0.0,5,6,8,0.039,0.17,0.179,6,4,8,0.017,0.092,0.079,7,5,7,0.032,0.161,0.153,8,4,7,0.01,0.085,0.088,9,1,4,0.0,0.058,0.0,潮流程序chaoliu.txt文档:#include<stdio.h>#include<math.h>#define N 9 /*总结点数*/#define M 6 /*PQ结点数*/#define K 9 /*线路数*/#define eps 1e-4void guass(int n,int m,float c[],float b[][N],float x[]) /*高斯函数*/{float a[N][N],y[N];int i,j,k;for(i=0;i<m;i++){y[i]=c[i];for(j=0;j<m;j++) a[i][j]=b[i+n][j+n];}for(k=0;k<m-1;k++)for(i=k+1;i<m;i++){for(j=k+1;j<m;j++) a[i][j]=a[i][j]-a[k][j]*a[i][k]/a[k][k]; y[i]=y[i]-y[k]*a[i][k]/a[k][k];}for(i=m-1;i>=0;i--){for(j=i+1;j<m;j++) y[i]-=a[i][j]*x[j];x[i]=y[i]/a[i][i];}}struct line{int Lindex;int Headnode;int Endnode;float R;float X;float b;}line[K];struct line *t;main(){float r[N]={0.0};float u[N]={1.04,1.025,1.025,1.0,1.0,1.0,1.0,1.0,1.0};float p[N]={1,1.63,0.85,0.0,0.0,0.0,-1.25,-0.9,-1.0};float q[N]={1,1,1,0.0,0.0,0.0,-0.5,-0.3,-0.35};float g[N][N]={0.0},b[N][N]={0.0};float h[N][N]={0.0};float B[N][N]={0.0};float temp;float H[K][6];float lr,lx,lb1,lg,lb;int i,j;int ku=0,kr=0,kp=1,kq=1;void val(float u[N],float g[N][N],float b[N][N],float r[N],int ku, int kr,float h[N][N]); /*函数申明*/FILE *fp;fp=fopen("data.txt","r");for(i=0;i<K;i++){for(j=0;j<6;j++){ fscanf(fp,"%f,",&temp);H[i][j]=temp;}}fclose(fp);for(i=0;i<K;i++){line[i].Lindex=(int)H[i][0];line[i].Headnode=(int)H[i][1] ;line[i].Endnode=(int)H[i][2];line[i].R=H[i][3];line[i].X=H[i][4];line[i].b=H[i][5];} for(t=line;t<line+K;t++){i=t->Headnode-1;j=t->Endnode-1;lr=t->R;lx=t->X;lb1=t->b;lg=lr/(lr*lr+lx*lx);lb=-lx/(lr*lr+lx*lx);g[i][i]+=lg;g[j][j]+=lg;b[i][i]+=lb+lb1;b[j][j]+=lb+lb1;h[i][j]=h[j][i]=-lb1;g[i][j]-=lg;g[j][i]-=lg;b[i][j]-=lb;b[j][i]-=lb;}getch();printf("\n=====jie dian dao na ju zhen=====\n");for(i=0;i<N;i++){for(j=0;j<N;j++){B[i][j]=b[i][j];printf("%8f,",B[i][j]);}printf("\n");}printf("\n");getch();printf("\n=====gei ding chu zhi=====\n");for(i=0;i<N;i++)printf("u[%d]=%8f p[%d]=%8f q[%d]=%8f\n",i+1, u[i],i+1,p[i],i+1,q[i]);printf("\n=====die dai qiu jie=====\n");while(kp==1){float ip,iq,max;float dp[N],dq[N],dr[N];float dpu[N],dqu[N];float y1[N-1],y2;float x1[N-1],x2;for(i=1;i<N;i++) /*算dp对应B',N-1维,除去平衡结点*/{ip=0;for(j=0;j<N;j++)ip=ip+u[j]*(g[i][j]*cos(r[i]-r[j])+b[i][j]*sin(r[i]-r[j ]));dp[i]=p[i]-u[i]*ip;dpu[i]=dp[i]/u[i];printf("dp[%d]=%8f\n",i+1,dp[i]);printf("\n");}getch();max=fabs(dpu[1]);for(i=1;i<N;i++){if(fabs(dpu[i])>max)max=fabs(dpu[i]);}if (max>=eps){for(i=0;i<N-1;i++)y1[i]=-dpu[i+1]; /*起值不同,为了对应,故加一*/guass(1,N-1,y1,B,x1);for (i=1;i<N;i++){dr[i]=x1[i-1]/u[i];r[i]=r[i]+dr[i];}printf("\n===== di %d ci die dai hou dian ya xiang jiao chu zhi =====\n",kr+1);for(i=1;i<N;i++)printf("r[%d]=%8f\n",i+1,r[i]);getch();printf("\n\n");kr=kr+1;kq=1;top:for(i=N-M;i<N;i++) /*算dq对应B",仅M维,除去平衡结点和PV 结点*/{iq=0;for(j=0;j<N;j++)iq=iq+u[j]*(g[i][j]*sin(r[i]-r[j])-b[i][j]*cos(r[i]-r [j]));dq[i]=q[i]-u[i]*iq;dqu[i]=dq[i]/u[i];printf("dq[%d]=%8f\n",i+1,dq[i]);printf("\n");}max=fabs(dqu[N-M]);for (i=N-M;i<N;i++){if(fabs(dqu[i])>max)max=fabs(dqu[i]);}if(max>=eps){for(i=0;i<M;i++)y2[i]=-dqu[i+N-M]; /*同上,对应加N-M*/guass(N-M,M,y2,B,x2);for(i=N-M;i<N;i++)u[i]=u[i]+x2[i-(N-M)];printf("\n=====di %d ci die dai dian ya chuzhi=====\n",ku+1);for(i=N-M;i<N;i++)printf("u[%d]=%8f\n",i+1,u[i]);printf("\n\n");ku=ku+1;kp=1;}else{kq=0;if(kp==0)val(u,g,b,r,ku,kr,h);}}else{kp=0;if(kq==0)val(u,g,b,r,ku,kr,h);elsegoto top;}}}void val(float u[N],float g[N][N],float b[N][N],float r[N],int ku, int kr,float h[N][N]){float ps=0,pv1=0,pv2=0;float qs=0,qv1=0,qv2=0;float p[N][N]={0};float q[N][N]={0};float s[N][N];float dp[N][N]={0};float dq[N][N]={0};float ds[N][N];float dSp=0,dSq=0;int i,j;FILE *fp1;printf("\n=====ping heng jie dian gong lv =====\n");getch();for(i=0;i<N;i++){ps=ps+u[0]*u[i]*(g[0][i]*cos(r[0]-r[i])+b[0][i]*sin(r[0 ]-r[i]));qs=qs+u[0]*u[i]*(g[0][i]*sin(r[0]-r[i])-b[0][i]*cos(r[0 ]-r[i]));}printf("sp=%8f+j(%8f)\n",ps,qs);printf("\n=====PV jie dian gong lv=====\n");getch();for(i=0;i<N;i++){pv1=pv1+u[1]*u[i]*(g[1][i]*cos(r[1]-r[i])+b[1][i]*sin(r [1]-r[i]));qv1=qv1+u[1]*u[i]*(g[1][i]*sin(r[1]-r[i])-b[1][i]*cos(r [1]-r[i]));}printf("sv1=%8f+j(%8f)\n",pv1,qv1);for(i=0;i<N;i++){pv2=pv2+u[2]*u[i]*(g[2][i]*cos(r[2]-r[i])+b[2][i]*sin(r [2]-r[i]));qv2=qv2+u[2]*u[i]*(g[2][i]*sin(r[2]-r[i])-b[2][i]*cos(r [2]-r[i]));}printf("sv2=%8f+j(%8f)\n",pv2,qv2);for(i=0;i<N;i++)for(j=0;j<N;j++){p[i][j]=u[i]*u[i]*(-g[i][j])+u[i]*u[j]*(g[i][j]*cos(r[ i]-r[j])+b[i][j]*sin(r[i]-r[j]));q[i][j]=u[i]*u[i]*(-h[i][j]+b[i][j])+u[i]*u[j]*(g[i][j ]*sin(r[i]-r[j])-b[i][j]*cos(r[i]-r[j]));dp[i][j]=u[i]*u[i]*(-g[i][j])+u[i]*u[j]*(g[i][j]*cos(r [i]-r[j])+b[i][j]*sin(r[i]-r[j]))+u[j]*u[j]*(-g[j][i])+u[j]*u[i]*(g[j][i]*cos(r[j]-r[i])+b[j][i]*sin(r[j]-r[i]));dq[i][j]=u[i]*u[i]*(-h[i][j]+b[i][j])+u[i]*u[j]*(g[i][ j]*sin(r[i]-r[j])-b[i][j]*cos(r[i]-r[j]))+u[j]*u[j]*(-h[j][i]+b[j][i])+u[j]*u[i]*(g[j][i]*sin(r[j]-r[i])-b[j][i]*cos(r[j]-r[i]));}printf("\n======die dai hou dian ya yu xiang jiao zhi======\n");getch();for(i=0;i<N;i++)printf("u[%d]=%8f r[%d]=%8f\n",i+1,u[i],i+1,r[i]); printf("\n=====xian lu gong lv=======\n");for(i=0;i<N;i++){getch();{ for(j=0;j<N;j++){ printf("s[%d][%d]=%8f+j(%8f)\n",i+1,j+1,p[i][j],q[i][j]); printf("\n");}printf("\n");}}printf("\n");printf("\n=====xian lu sun hao=====\n");for(i=0;i<N;i++){getch();{ for(j=i+1;j<N;j++){ printf("ds[%d][%d]=%8f+j%8f\n",i+1,j+1,dp[i][j],dq[i][j]); printf("\n");}printf("\n");}}printf("\n");printf("\n=====wang luo zong sun hao=====\n");getch();for(i=0;i<N;i++){for(j=i+1;j<N;j++){ dSp+=dp[i][j];dSq+=dq[i][j];}}printf("dS=%8f+j(%8f)\n",dSp,dSq);printf("\n======die dai ci shu========\n");printf("ku=%d\n",ku);printf("kr=%d\n",kr);printf("\n=======shu ju bao cun=====\n");fp1=fopen("jieguo.txt","w+");{fprintf(fp1,"xian lu cao liu:\n");for(i=0;i<N;i++){for(j=0;j<N;j++)fprintf(fp1,"s[%d][%d]=%8f+j%8f\n",i+1,j+1,p[i][j],q[i][j]); }fprintf(fp1,"\n");fprintf(fp1,"zhi lu sun hao:\n");for(i=0;i<N;i++){for(j=i+1;j<N;j++)fprintf(fp1,"ds[%d][%d]=%8f+j%8f\n",i+1,j+1,dp[i][j],dq[i][j ]);}fprintf(fp1,"\n");fprintf(fp1,"wang luo zong sun hao:\n");fprintf(fp1,"dS=%8f+j(%8f)\n",dSp,dSq);fprintf(fp1,"\n");fprintf(fp1,"PV jie dian gong lv:\n");fprintf(fp1,"sv1=%8f+j(%8f)\n",pv1,qv1);fprintf(fp1,"sv2=%8f+j(%8f)\n",pv2,qv2);fprintf(fp1,"\n");fprintf(fp1,"ping heng jie dian gong lv:\n");fprintf(fp1,"sp=%8f+j(%8f)\n",ps,qs);fprintf(fp1,"\n");fprintf(fp1,"jie dian dian ya yu xiang jiao:\n");for(i=0;i<N;i++)fprintf(fp1,"u[%d]=%8f r[%d]=%8f\n",i+1,u[i],i+1, r[i]);}printf("\n===========THE END==============\n");getch();}。