电力系统潮流计算C语言程序及说明

实验目的根据所给的电力系统，编制潮流计算程序，通过计算机进行调试，最后完成一个切实可行的电力系统计算应用程序。通过自己设计电力系统计算程序使同学对电力系统分析有进一步理解，同时加强计算机实际应用能力的训练。程序计算原理1、概述应用计算机进行电力系统计算，首先要掌握电力系统相应计算的数学模型；其次是运用合理的计算方法；第三则是选择合适的计算机语言编制计算程序。建立电力系统计算的相关数学模型，就是建立用于描述电力系统相应计算的有关参数间的相互关系的数学方程式。该数学模型的建立往往要突出问题的主要方面，即考虑影响问题的主要因素，而忽略一些次要因素，使数学模型既能正确地反映实际问题，又使计算不过于复杂。运用合理的计算方法，就是要求所选用的计算方法能快速准确地得出正确结果，同时还应要求在解算过程中占用内存少，以利提高计算机的解题规模。选择合适的语言编写程序，就是首先确定用什么计算机语言来编制程序；其次是作出计算的流程图；第三根据流程图用选择的语言编写计算程序。然后上机调试，直到语法上无错误。本程序采用C语言进行编程。所编制的程序难免存在逻辑错误，因此先用一个已知结果的系统作为例题进行计算。用程序计算的结果和已知结果相比较，如果结果相差甚远就要逐步分析程序的计算步骤，查出问题的出处；如果结果比较接近，则逐步分析误差来源；直到结果正确为止。2、电力系统潮流计算的程序算法潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压（幅值及相角）、网络中的功率分布及功率损耗等。目前计算机潮流计算的方法主要有牛顿-拉夫逊算法和PQ分解法。牛顿-拉夫逊算法是数学上求解非线形方程组的有效方法，具有较好的收敛性，曾经是潮流计算中应用比较普遍的方法。PQ快速分解法是从牛顿-拉夫逊算法演变而来的，是将纯数学的牛顿-拉夫逊算法与电力系统具体特点相结合并进行简化与改进而得出的。PQ快速分解法比牛顿-拉夫逊算法大大提高了计算速度和节省了内存，故而本程序以PQ快速分解法进行潮流计算。1）形成节点导纳矩阵（1）自导纳的形成对节点i其自导纳Yii是节点i以外的所有节点都接地时节点i对地的总导纳。显然，Yii应等于与节点i相接的各支路导纳之和，即Y二y+Xyiii0 ijj式中，yi0为节点i与零电位节点之间的支路导纳；yij为节点i与节点j之间的支路导纳。2）互导纳的形成对节点i与节点k之间的互导纳是节点i、k之间的支路导纳的负值，即y二-yikik不难理解Y二Y。若节点i和k没有支路直接相连时，便有Yik=°kiik（3）含变压器支路的处理若节点p、q间接有变压器，如下图所示，则可作出其n型等值电路为：图1变压器口型等值电路则p、q的自导纳和节点间的互导纳分别为2） 计算不平衡功率AP'AQ并形成修正方程式对每一个PQ节点或每一个PV节点都可以根据下列公式计算出有功功率增量AP而对于每一个PQ节点还可以根据下面的公式计算出无功功率增量AQ在有功功率增量和无功功率增量不满足如下约束条件时利用PQ分解法则可以形成如下修正方程3） 利用因子表法求解修正方程在电网计算中经常遇到这样的问题，对方程组需要反复多次求解，而每次求解仅改变常数项F,系数矩阵保持不变。按照一般的高斯消去法，对每一改变的常数项，形成包括常数项及系数矩阵在内的增广矩阵，然后消去回代求出其解。可以看出，每次对增广矩阵中A矩阵元素的消元都是重复的，为了避免这种重复，我们把对相同的系数矩阵重复进行的消去与对不同的常数项进行的消去分开进行，因此对系数矩阵的消去只需进页脚内容2行一次，并在消去的过程中将对常数项进行消去运算的运算因子保存下来，形成所谓因子表，这就是因子表法。因为因子表记录了高斯消去法对常数项进行消去的全部信息，利用它便可对不同常数项进行消去，形成上三角矩阵，最后求出全部未知数。在使用PQ分解法时，其系数矩阵是在迭代过程中保持不变的，所以为了节省内存和缩短运算时间我们采取了因子表法。同时由于电网的节点导纳矩阵矩阵是稀疏阵和对称阵，于是我们可以采取只保存系数矩阵的上三角阵来使运算更为简化。aaaLa111213L1naaa22232naLa333nOM其中A=nn称为系数矩阵，其中A=nn称为系数矩阵，aX=[xxxLxA称为未知数向量，123nF=[fffLf｝称为常数项向量。将矩阵A的元素进行如下处理:123n得到因子表其中D=1a(i-i),U=a(i)(i<j)；ii ii ijij再利用因子表进行前代过程，求出每次迭代后的常数项。其前代公式是求得向量F=rf⑴f⑵f⑶Lf(n)T；1 2 3 n再由因子表与前代得到的向量F，得到方程组求解出此方程即可得到线性方程组的解向量X=[xxxLxb。1 2 3n4）多次迭代最终求得V和6以及全线路功率利用上面所介绍的方法求解修正方程组可以求得VV二[vvVVVVLVV]T和V6=[v6V6V6LV6]T。12 3 m 1 23 n-1再利用6（k+1）=6（k）+V6仏），V（k+1）=V（k）+VV（k）求得每次迭代后的结果。多次迭代当其满足iiiiii约束条件max｛P（k）｝＜e和max[vQ（k）｝＜e时，迭代结束。迭代结束后即可得到各节点的Vi P i Q和6，再根据V、6来计算PV节点的无功功率Q和平衡节点的功率以及网络中的功率分布。PV节点及平衡节点无功功率计算公式为：平衡节点有功功率计算公式为：以下图所标示的正方向，输电线路功率的计算公式如下图2支路功率计算对其进行实部虚部进行分解可得P、Q计算公式为：程序及说明1、主要变量说明1）结构体类型说明1）节点功率结构体structNodetype{floatP,Q}其中，P为节点的有功功率，Q为无功功率。节点功率不区分负荷功率和发电机功率，其值为本节点连接的各支路输入功率及节点所接负荷、发电机功率之和，且规定功率流入节点为正，流出为负。详细说明参见下一章“算例及结果”的第二节“源数据格式说明”。（2）线路参数结构体structLinetype{floatG,B,B0,k}其中，G、B为线路的导纳和容纳；B0为线路的考虑变压器n型等值电路后的对地充电容纳的一半Bc/2；k为折算到标准变压器支路后的变压器变比。详细说明参见下一章“算例及结果”的第二节“源数据格式说明”。2）变量说明

表2程序主要变量说明主要变量类型含义Nodeint系统总节点数NPintPV+PQ节点数，即非平衡节点数NQintPQ节点数Numint*原始节点编号与程序表示编号映射数组NostructNodetype*节点功率数组Vfloat*节点电压有效值数组Ditafloat*节点电压相角数组YstructLinetype**线路参数矩阵BP、BQfloat**有功、无功简化雅克比矩阵B、B”countunsignedintPQ迭代次数eP、eQconst有功、无功迭代精度控制kp、kqint有功、无功迭代结束标志dP、dQfloat*有功、无功不平衡量数组3、程序流程图图3程序主流程图图4迭代部分流程图4、程序代码/*FUNCTION:POWERFLOW*//*WRITTENBY:HUANG&YANG&TONG*//*LASTEDITED:2008-11-24*/#include<stdio.h>#include<math.h>/***宏定义***/#defineeP0.00001#defineeQ0.00001#defineY_(i,j)(*(*(Y+i)-i+j))#defineYij(*(Yi+j))#defineYji(*(*(Y+j)-j+i))#definePjiYji.G*cos(tmp)+Yji.B*sin(tmp)#definePijYij.G*cos(tmp)+Yij.B*sin(tmp)#defineQjiYji.G*sin(tmp)-Yji.B*cos(tmp)#defineQijYij.G*sin(tmp)-Yij.B*cos(tmp)/***结构体变量定义***/structNodetype/*节点功率*/{floatP,Q}structLinetype/*线路类型*/{floatG,B,B0,k}/***子函数声明***/voidin_node()/*读节点功率*/voidin_line()/*读线路参数*/voidB_Form()/*生成BP、BQ矩阵*/voidfactor()/*求因子表*/voidsolve(float**B,float*X,intN)/*解方程组*/voidPrtNode()/*打印节点参数*/voidErrorMsg(intFlag)/*错误提示信息*//***全局变量声明***/intNode/*节点数*/int*Num/*保存原始节点序号*/intNP,NQ=0/*PV+PQ、PQ节点数*/structNodetype*No/*节点数据*/structLinetype**Y/*线路参数*/float**BP,**BQ/*有功、无功简化雅克比矩阵B*/float*V/*节点电压有效值*/float*Dlta/*节点电压相角值*/unsignedintcount=0/*迭代计数*/inti,j,k,m/*通用下标值*/floattmp/*临时数据暂存*/char*Type/*节点类型*/FILE*in,*out/*输入、输出文件指针*//**I****I****I**主函数**I****I****I**/intmain(void){intkp=1,kq=1/*P、Q精度判断：1-不满足，0-满足精度要求*/float*dP,*dPi,*dQ,*dQi/*AP、AQ*/floatDltaistructLinetype*YistructNodetype*NoifloattP,tQif((in=fopen("Data.txt","r"))==NULL)ErrorMsg(1)if((out=fopen("out.txt","w"))==NULL)ErrorMsg(2)in_node()/*读取节点参数并统计PQ、PV节点数*/in_line()/*读取线路参数并形成Y矩阵*/B_Form()/*形成B(BP&BQ)矩阵*/factor。/*求B因子式(仍保存于BP&BQ)*/for(i=Oi<NQi++)*(V+i)=l/*对PQ节点电压V赋初值*/dP=(float*)malloc(sizeof(float)*NP)/*AP*/dQ=dP/*AQ*//*AP、AQ不同时存在，故而可共用空间*/loop：/****开始迭代****/if(kp==O&&kq==O)gotoloopEndcount++/*迭代次数加一*/if(count==65535)ErrorMsg(99)/*不收敛，退出*/kp=Ofor(i=Oi<NPi++){dPi=dP+iYi=*(Y+i)-iDltai=*(Dlta+i)*dPi=0for(j=0j<Nodej++){tmp=Dltai-*(Dlta+j)/*tmp即§ij*/if(i>j)*dPi+=*(V+j)*(Pji)else*dPi+=*(V+j)*(Pij)}/*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/*dPi*=*(V+i)*dPi=(*(No+i)).P-*dPi/*求得APi*/if(fabs(*dPi)>0x8fffffff)ErrorMsg(99)/*不收敛，退出*/if(fabs(*dPi)>eP)kp=l/*有不满足精度的AP即令kp=1*/*dPi/=*(V+i)/*求得常数项APi/Vi*/}if(kp==0)gotoloopQsolve(BP,dP,NP)for(i=Oi<NPi++)*(Dlta+i)+=(*(dP+i)/(*(V+i)))/*修正相角6+=AS*/loopQ:if(kp==O&&kq==O)gotoloopEndkq=Ofor(i=Oi<NQi++){dQi=dQ+iYi=*(Y+i)-iDltai=*(Dlta+i)*dQi=Ofor(j=Oj<Nodej++){tmp=Dltai-*(Dlta+j)/*tmp即6ij*/if(i>j)*dQi+=*(V+j)*(Qji)else*dQi+=*(V+j)*(Qij)}/*注意到y矩阵为上三角矩阵，i>j时要交换下标*/*dQi*=*(V+i)*dQi=(*(No+i)).Q-*dQi/*求得AQi*/if(fabs(*dQi)>0x8fffffff)ErrorMsg(99)/*不收敛，退出*/f(fabs(*dQi)>eQ)kq=l/*有不满足精度的AQ即令kq=l*/*dQi/=*(V+i)/*求得常数项AQi/Vi*/}if(kq==0)gotoloopsolve(BQ,dQ,NQ)for(i=0i<NQi++)*(V+i)+=*(dQ+i)/*修正PQ节点电压V+=AV*/gotoloop/*无功迭代，则必定需要下一轮回迭代判断*/loopEnd：/****迭代结束****/free(dP)/*释放内存空间*//****计算PV节点和平衡节点的无功功率Q****/for(i=NQi<Nodei++){Noi=No+iYi=*(Y+i)-iDltai=*(Dlta+i)for(j=0j<Nodej++){tmp=Dltai-*(Dlta+j)/*tmp即§ij*/if(i>j)(*Noi).Q+=*(V+j)*(Qji)else(*Noi).Q+=*(V+j)*(Qij)}/*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/(*Noi).Q*=*(V+i)}/****计算平衡节点的有功功率P****/i=NPNoi=No+iDltai=*(Dlta+i)for(j=0j<Nodej++){tmp=Dltai-*(Dlta+j)/*tmp即§ij*/(*Noi).P+=*(V+j)*(Pji)}/*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/(*Noi).P*=*(V+i)/****输出最终结果****/fprintf(out,"\n\n【潮流计算结果(节点)】(迭代次数k=%3d)\n",count-l)PrtNode()/****计算全部线路功率****/fprintf(out,"\n\n【潮流计算结果(线路)】\n")fprintf(out,"线路PQ\n")for(k=0k<Nodek++){i=*(Num+k)Yi=*(Y+i)-iDltai=*(Dlta+i)Noi=No+ifor(m=0m<Nodem++){j=*(Num+m)if(j==i)continuetmp=Dltai-*(Dlta+j)/*tmp即§ij*/if(j<i){f(Yji.B==0)continue/*若Bij=O,则节点i、j无直接联系*/tP=*(V+j)*(Pji)tP=*(V+i)*Yji.G-tPtP*=*(V+i)tQ=-*(V+j)*(Qji)tQ-=*(V+i)*(Yji.B-Yji.B0/Yji.k)tQ*=*(V+i)}else{f(Yij.B==0)continue/*若Bij=O,则节点i、j无直接联系*/tP=*(V+j)*(Pij)tP=*(V+i)*Yij.G-tPtP*=*(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)}}fclose(out)system("cmd/cstartout.txt")return(0)/**f****f****f**主函数**f****f****f**//******************主函数******************//****************子函数：读节点数据****************/voidin_node(){structNodetype*Noi/*临时变量*/fscanf(in,"%d%d",&Node,&k)/*读取节点数Node*/NP=Node-1/*PV+PQ节点数，即非平衡节点数目*/Num=(int*)malloc(sizeof(int)*Node)/*开Node个空间，每节点一个*/V=(float*)malloc(sizeof(float)*Node)/*电压*/Dlta=(float*)malloc(sizeof(float)*Node)/*电压相角*/No=(structNodetype*)malloc(sizeof(structNodetype)*Node)/*节点功率*/j=1while(k!=O)/*若k=0,表明节点数据读取完毕*/{switch(k){case1:k=NQNQ++break/*NQ统计PQ节点个数*/case2:k=NP-jj++break/*从NP-1个空间倒着保存PV节点*/case3:k=NPbreak/*平衡节点*/default:ErrorMsg(3)}Noi=No+kfscanf(in,"%d%f%f%f%f",&i,&(*Noi).P,&(*Noi).Q,V+k,Dlta+k)i--/*节点编号减一，以和程序表达方式兼容*/*(Num+i)=k/*第i个Num元素中存放i节点在No中的下标*/fscanf(in,"%d",&k)/*读取节点类型*/}if(NQ+j!=Node)ErrorMsg(4)/*检验节点数据个数是否够Node个*/fprintf(out,"【节点参数表】\n")PrtNode()fprintf(out,"总节点：％d\nPQ节点：％d\nPV节点：％d\n",Node,NQ,NP-NQ)}/************子函数：读线路数据，并形成节点导纳矩阵Y************/voidin_line(){structLinetype*YifloatR,X,k,Bm=sizeof(structLinetype)Y=(structLinetype**)malloc(m*Node)/*先开Node行,每一个节点一行*/for(i=0i<Nodei++)/*再在第i行上面开辟Node-i列*/{/*即以上三角存储Y矩阵*/*(Y+i)=(structLinetype*)malloc(m*(Node-i))Yi=*(Y+i)-ifor(j=ij<Nodej++){Yij.G=Yij.B=Yij.B0=Yij.k=0}/*初始化*/}while(feof(in)==0)/*文件指针到文件末*/{fscanf(in,"%d%d%f%f%f%f",&i,&j,&R,&X,&k,&B)i--j--i=*(Num+i)/*转换节点号为该节点在程序中的储存编号*/j=*(Num+j)(*(*(Y+i))).B+=B/*将对地充电导纳累加到自导纳*/(*(*(Y+j))).B+=Bif(k!=1.0){X*=kR=0tmp=(1-k)/X(*(*(Y+i))).B+=tmp/*将变压器的对地充电容纳累加到自导纳*/(*(*(Y+j))).B+=-(tmp/k)B=tmpk=-k}if(i>j){tmp=ii=jj=tmpk=1/kB*=k}Yi=*(y+i)-i/*以Yi代替*(Y+i)-i,简化表达式并避免重复计算*/Yij.BO二B/*保存ijO、jiO对地充电电容到BijO*/Yij.k=k/*且有B0ji=B0ij/k*/tmp=R*R+X*XR/=tmpX/=tmpYij.G=-R/*生成互导纳*/Yij.B=X(*(*(Y+i))).G+=R/*将线路互导纳累加到自导纳*/(*(*(Y+i))).B+=-X(*(*(Y+j))).G+=R(*(*(Y+j))).B+=-X}fclose(in)fprintf(out,"\n【节点导纳矩阵Y】\n")for(k=0k<Nodek++){i=ki=*(Num+i)/*查取第i节点在程序中存储序号*/for(j=0j<kj++)fprintf(out,"\t\t\t")for(m=km<Nodem++){j=*(Num+m)/*查取第j节点在程序中存储序号*/if(i<j)fprintf(out,"(%10.6f,%10.6f)",Y_(i,j).G,Y_(i,j).B)elsefprintf(out,"(%10.6f,%10.6f)",Y_(j,i).G,Y_(j,i).B)}fprintf(out,"\n")}/****************子函数：生成BP、BQ矩阵****************/voidB_Form()float*BPi,*BQistructLinetype*Yiintsize=sizeof(float)BP=(float**)malloc(size*NP)/*以上三角存储*/for(i=0i<NPi++)*(BP+i)=(float*)malloc(size*(NP-i))for(i=0i<NPi++){BPi=*(BP+i)-i/*以BPi代替*(BP+i)-i,避免重复计算*/Yi=*(Y+i)-ifor(j=ij<NPj++)*(BPi+j)=Yij.B/*(BPi+j)即相当于BP[i][j]*/}BQ=BP/*BP包含BQ,BP左上角的NQ*NQ子阵即BQ*//****************子函数：求因子表****************/voidfactor(){float*BPi,*BPk,*BQifor(i=0i<NPi++){BPi=*(BP+i)-ifor(k=0k<ik++){BPk=*(BP+k)-ktmp=(*(BPk+i))/(*(BPk+k))for(j=ij<NPj++)(*(BPi+j))-=tmp*(*(BPk+j))}*(BPi+i)=1/(*(BPi+i))for(j=i+1j<NPj++)*(BPi+j)*=*(BPi+i)}/****************子函数：解方程组****************/voidsolve(float**B,float*X,intN){float*Bi,*Xifor(i=0i<Ni++)*(X+i)=-*(X+i)/*对常数项取负*//****对常数列进行前代****/for(i=0i<Ni++){Bi=*(B+i)-iXi=X+ifor(j=i+1j<Nj++)*(X+j)-=*(Bi+j)**Xi*Xi*=*(Bi+i)}/****回代以求解方程组****/for(i=N-1i>=0i--){Bi=*(B+i)-iXi=X+ifor(j=N-1j>ij--)*Xi-=*(Bi+j)**(X+j)}}/****打印节点参数****/voidPrtNode(){structNodetype*Noifprintf(out,"节点类型PQV6\n")for(i=0i<Node){j=*(Num+i)/*查取第i节点在程序中存储序号*/Noi=No+jif(j<NQ)Type="PQ"elseType="PV"if(j==NP)Type="BS"fprintf(out,"%3d%s%10.6f%10.6f%10.6f%10.6f\n",++i,Type,(*Noi).P,(*Noi).Q,*(V+j),*(Dlta+j)/0.017453)}/****************子函数：错误信息****************/voidErrorMsg(intFlag){switch(Flag){case1:printf("\n\tError(1):FailedtoOpenFile\"Data.txt\"!")breakcase2:printf("\n\tError(2):FailedtoCreatFile\"out.txt\"!")breakcase3:printf("\n\tError(3):NodeDataError,PleaseCheck!")breakcase4:printf("\n\tError(4):LackNodeData,PleaseCheck!")breakcase99:printf("k=%d\n\tError(99):It'sEmanative!",count)break}getch()fclose(out)exit(Flag)}算例及结果1、原始数据程序编写过程中，采用了《电力系统分析》教材的4节点、IEEE9节点参数进行了验证。本说明书中仅列出4节点验算过程。4节点测试数据来源于参考文献中的《电力系统分析》(下册)第61页例11-5，用于检验程序正确性的参考计算结果见67页。对其基本情况总结如下。网络连接图图6算例电力系统网络图节点参数表4算例节点参数表

节点类型节点编号有功功率P无功功率Q电压V电压相角811-0.30-0.180012-0.55-0.1300230.501.10034001.0503）线路参数表5算例线路参数表节点1节点2阻抗R感抗X变压器变比k充电容纳Bc/2120.10.410.01528140.120.510.01920240.080.410.014133100.30.90909104）数据文件data.txt434001.050230.501.10011-0.30-0.1800图7数据文件格式2、源数据格式说明源数据有功功率、无功功率、电压、阻抗、感抗、对地充电容纳均以标幺值表示。数据文件必须命名为data.txt且与潮流计算程序放置于同一个文件中。数据文件data.txt包含两类参数：节点参数和线路参数。节点数据块与线路数据块之间用数字0作为间隔，即在节点数据块结束后，另起一行输入0，然后再在后面按格式要求录入线路参数1）节点参数节点参数包括：系统节点数Node、节点功率（有功P、无功Q）,节点电压（有效值页脚内容35V、相角Delta）参数组织格式：（1）节点数Node节点数Node写在参数文件的开头，如：4表明为四节点系统。（2） 功率和电压P/Q/V/Delta首先给出节点参数示例：230.501.100第一列数字2表明该节点的类型为2-PV节点；第二列数字3表明该行数据为节点3的参数；后面三列依次为P、Q、V的给定值，给定值为0，表明该项参数未知；第六列为相角非平衡节点的§即为PQ迭代的初始相角值，平衡节点的即为给定的相角值。节点功率为各支路输入功率之和，且规定功率流入节点为正，流出为负。故而负荷功率为负值，发电机功率为正值。（3） 节