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

1、电力系统潮流计算C语言程序及说明实验目的根据所给的电力系统，编制潮流计算程序，通过计算机进行调试，最后完成一个切实可行的电力系统计算应用程序。通过自己设计电力系统计算程序使同学对电力系统分析有进一步理解，同时加强计算机实际应用能力的训练。程序计算原理1、概述应用计算机进行电力系统计算，首先要掌握电力系统相应计算的数学模型；其次是运用合理的计算方法；第三则是选择合适的计算机语言编制计算程序。建立电力系统计算的相关数学模型，就是建立用于描述电力系统相应计算的有关参数间的相互关系的数学方程式。该数学模型的建立往往要突出问题的主要方面，即考虑影响问题的主要因素，而忽略一些次要因素，使数学模型既能正确地

2、反映实际问题，又使计算不过于复杂。运用合理的计算方法，就是要求所选用的计算方法能快速准确地得出正确结果，同时还应要求在解算过程中占用内存少，以利提高计算机的解题规模。选择合适的语言编写程序，就是首先确定用什么计算机语言来编制程序；其次是作出计算的流程图；第三根据流程图用选择的语言编写计算程序。然后上机调试，直到语法上无错误。本程序采用C语言进行编程。所编制的程序难免存在逻辑错误，因此先用一个已知结果的系统作为例题进行计算。用程序计算的结果和已知结果相比较，如果结果相差甚远就要逐步分析程序的计算步骤，查出问题的出处；如果结果比较接近，则逐步分析误差来源；直到结果正确为止。2、电力系统潮流计算的程

3、序算法潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压（幅值及相角）、网络中的功率分布及功率损耗等。目前计算机潮流计算的方法主要有牛顿-拉夫逊算法和PQ分解法。牛顿-拉夫逊算法是数学上求解非线形方程组的有效方法，具有较好的收敛性，曾经是潮流计算中应用比较普遍的方法。PQ快速分解法是从牛顿-拉夫逊算法演变而来的，是将纯数学的牛顿-拉夫逊算法与电力系统具体特点相结合并进行简化与改进而得出的。PQ快速分解法比牛顿-拉夫逊算法大大提高了计算速度和节省了内存，故而本程序以PQ快速分解法进行潮流计算。1）形成节点导纳矩阵（1）自导纳的形成对节点i其自

4、导纳Yii是节点i以外的所有节点都接地时节点i对地的总导纳。显然，Yii应等于与节点i相接的各支路导纳之和，即式中，yi0为节点i与零电位节点之间的支路导纳；yij为节点i与节点j之间的支路导纳。（2）互导纳的形成对节点i与节点k之间的互导纳是节点i、k之间的支路导纳的负值，即不难理解。若节点i和k没有支路直接相连时，便有Yik=0（3）含变压器支路的处理若节点p、q间接有变压器，如下图所示，则可作出其型等值电路为：图1 变压器型等值电路则p、q的自导纳和节点间的互导纳分别为2）计算不平衡功率P、Q并形成修正方程式对每一个PQ节点或每一个PV节点都可以根据下列公式计算出有功功率增量P而对于每一

5、个PQ节点还可以根据下面的公式计算出无功功率增量Q在有功功率增量和无功功率增量不满足如下约束条件时利用PQ分解法则可以形成如下修正方程3）利用因子表法求解修正方程在电网计算中经常遇到这样的问题，对方程组需要反复多次求解，而每次求解仅改变常数项F，系数矩阵保持不变。按照一般的高斯消去法，对每一改变的常数项，形成包括常数项及系数矩阵在内的增广矩阵，然后消去回代求出其解。可以看出，每次对增广矩阵中A矩阵元素的消元都是重复的，为了避免这种重复，我们把对相同的系数矩阵重复进行的消去与对不同的常数项进行的消去分开进行，因此对系数矩阵的消去只需进行一次，并在消去的过程中将对常数项进行消去运算的运算因子保存下

6、来，形成所谓因子表，这就是因子表法。因为因子表记录了高斯消去法对常数项进行消去的全部信息，利用它便可对不同常数项进行消去，形成上三角矩阵，最后求出全部未知数。在使用PQ分解法时，其系数矩阵是在迭代过程中保持不变的，所以为了节省内存和缩短运算时间我们采取了因子表法。同时由于电网的节点导纳矩阵是稀疏阵和对称阵，于是我们可以采取只保存系数矩阵的上三角阵来使运算更为简化。若线性方程组一般形式如下：其中称为系数矩阵，称为未知数向量， 称为常数项向量。将矩阵A的元素进行如下处理：得到因子表其中 ；再利用因子表进行前代过程，求出每次迭代后的常数项。其前代公式是：求得向量；再由因子表与前代得到的向量F，得到方

7、程组求解出此方程即可得到线性方程组的解向量。4）多次迭代最终求得V和以及全线路功率利用上面所介绍的方法求解修正方程组可以求得和。再利用求得每次迭代后的结果。多次迭代当其满足约束条件和时，迭代结束。迭代结束后即可得到各节点的V和，再根据V、来计算PV节点的无功功率Q和平衡节点的功率以及网络中的功率分布。PV节点及平衡节点无功功率计算公式为：平衡节点有功功率计算公式为：以下图所标示的正方向，输电线路功率的计算公式如下：图2 支路功率计算对其进行实部虚部进行分解可得P、Q计算公式为：程序及说明1、主要变量说明1）结构体类型说明（1）节点功率结构体struct Nodetype float P,Q;其

8、中，P为节点的有功功率，Q为无功功率。节点功率不区分负荷功率和发电机功率，其值为本节点连接的各支路输入功率及节点所接负荷、发电机功率之和，且规定功率流入节点为正，流出为负。详细说明参见下一章“算例及结果”的第二节“源数据格式说明”。（2）线路参数结构体struct Linetype float G,B,B0,k;其中，G、B为线路的导纳和容纳；B0为线路的考虑变压器型等值电路后的对地充电容纳的一半Bc/2；k为折算到标准变压器支路后的变压器变比。详细说明参见下一章“算例及结果”的第二节“源数据格式说明”。2）变量说明表2 程序主要变量说明主要变量类型含义Nodeint系统总节点数NPintPV

9、+PQ节点数，即非平衡节点数NQintPQ节点数Numint *原始节点编号与程序表示编号映射数组Nostruct Nodetype *节点功率数组Vfloat *节点电压有效值数组Dltafloat *节点电压相角数组Ystruct Linetype *线路参数矩阵BP、BQfloat *有功、无功简化雅克比矩阵B、B”countunsigned intPQ迭代次数eP、eQconst有功、无功迭代精度控制kp、kqint有功、无功迭代结束标志dP、dQfloat *有功、无功不平衡量数组3、程序流程图图3 程序主流程图图4 迭代部分流程图4、程序代码/* FUNCTION : POWER

10、FLOW */* WRITTEN BY : HUANG&YANG&TONG */* LAST EDITED : 2008-11-24 */#include #include /* 宏定义 */#define eP 0.00001#define eQ 0.00001#define Y_(i,j) (*(*(Y+i)-i+j)#define Yij (*(Yi+j)#define Yji (*(*(Y+j)-j+i)#define Pji Yji.G*cos(tmp)+Yji.B*sin(tmp)#define Pij Yij.G*cos(tmp)+Yij.B*sin(tmp)#define Qj

11、i Yji.G*sin(tmp)-Yji.B*cos(tmp)#define Qij Yij.G*sin(tmp)-Yij.B*cos(tmp)/* 结构体变量定义 */struct Nodetype /* 节点功率 */ float P,Q;struct Linetype /* 线路类型 */ float G,B,B0,k;/* 子函数声明 */void in_node(); /* 读节点功率 */void in_line(); /* 读线路参数 */void B_Form(); /* 生成BP、BQ矩阵 */void factor(); /* 求因子表 */void solve(float

12、 *B,float *X,int N); /* 解方程组 */void PrtNode(); /*打印节点参数*/void ErrorMsg(int Flag); /* 错误提示信息 */* 全局变量声明 */int Node; /*节点数*/int *Num; /*保存原始节点序号*/int NP,NQ=0; /*PV+PQ、PQ节点数*/struct Nodetype *No; /*节点数据*/struct Linetype *Y; /*线路参数*/float *BP,*BQ; /*有功、无功简化雅克比矩阵B*/float *V; /*节点电压有效值*/float *Dlta; /*节点电

13、压相角值*/unsigned int count=0; /*迭代计数*/int i,j,k,m; /*通用下标值*/float tmp; /*临时数据暂存*/char *Type; /*节点类型*/FILE *in,*out; /*输入、输出文件指针*/* 主函数 */* 主函数 */int main(void) int kp=1,kq=1; /*P、Q精度判断：1-不满足，0-满足精度要求*/ float *dP,*dPi,*dQ,*dQi; /*P、Q*/ float Dltai; struct Linetype *Yi; struct Nodetype *Noi; float tP,tQ

14、; if(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=0;iNQ;i+) *(V+i)=1; /*对PQ节点电压V赋初值*/ dP=(float *)malloc(sizeof(float)*NP); /*P*/ dQ=dP;

15、/*Q*/*P、Q不同时存在，故而可共用空间*/loop: /*开始迭代*/ if(kp=0&kq=0) goto loopEnd; count+; /*迭代次数加一*/ if(count=65535) ErrorMsg(99); /*不收敛，退出*/ kp=0; for(i=0;iNP;i+) dPi=dP+i; Yi=*(Y+i)-i; Dltai=*(Dlta+i); *dPi=0; for(j=0;jj) *dPi+=*(V+j)*(Pji); else *dPi+=*(V+j)*(Pij); /*注意到Y矩阵为上三角矩阵，ij时要交换下标*/ *dPi*=*(V+i); *dPi=(

16、*(No+i).P-*dPi; /*求得Pi*/ if(fabs(*dPi)0x8fffffff) ErrorMsg(99); /*不收敛，退出*/ if(fabs(*dPi)eP) kp=1; /*有不满足精度的P即令kp=1*/ *dPi/=*(V+i); /*求得常数项Pi/Vi*/ if(kp=0) goto loopQ; solve(BP,dP,NP); for(i=0;iNP;i+) *(Dlta+i)+=(*(dP+i)/(*(V+i); /*修正相角+=*/loopQ: if(kp=0&kq=0) goto loopEnd; kq=0; for(i=0;iNQ;i+) dQi=

17、dQ+i; Yi=*(Y+i)-i; Dltai=*(Dlta+i); *dQi=0; for(j=0;jj) *dQi+=*(V+j)*(Qji); else *dQi+=*(V+j)*(Qij); /*注意到Y矩阵为上三角矩阵，ij时要交换下标*/ *dQi*=*(V+i); *dQi=(*(No+i).Q-*dQi; /*求得Qi*/ if(fabs(*dQi)0x8fffffff) ErrorMsg(99); /*不收敛，退出*/ if(fabs(*dQi)eQ) kq=1; /*有不满足精度的Q即令kq=1*/ *dQi/=*(V+i); /*求得常数项Qi/Vi*/ if(kq=0

18、) goto loop; solve(BQ,dQ,NQ); for(i=0;iNQ;i+) *(V+i)+=*(dQ+i); /*修正PQ节点电压V+=V*/ goto loop; /*无功迭代，则必定需要下一轮回迭代判断*/loopEnd: /*迭代结束*/free(dP); /*释放内存空间*/*计算PV节点和平衡节点的无功功率Q*/ for(i=NQ;iNode;i+) Noi=No+i; Yi=*(Y+i)-i; Dltai=*(Dlta+i); for(j=0;jj) (*Noi).Q+=*(V+j)*(Qji); else (*Noi).Q+=*(V+j)*(Qij); /*注意到

19、Y矩阵为上三角矩阵，ij时要交换下标*/ (*Noi).Q*=*(V+i); /*计算平衡节点的有功功率P*/ i=NP; Noi=No+i; Dltai=*(Dlta+i); for(j=0;jj时要交换下标*/ (*Noi).P*=*(V+i);/*输出最终结果*/fprintf(out,nn【 潮流计算结果(节点) 】 ( 迭代次数k=%3d )n,count-1);PrtNode();/*计算全部线路功率*/fprintf(out,nn【 潮流计算结果(线路) 】n);fprintf(out, 线路 P Qn); for(k=0;kNode;k+) i=*(Num+k); Yi=*(Y

20、+i)-i; Dltai=*(Dlta+i); Noi=No+i; for(m=0;mNode;m+) j=*(Num+m); if(j=i) continue; tmp=Dltai-*(Dlta+j); /*tmp即ij*/ if(ji) if(Yji.B=0) continue; /*若Bij=0，则节点i、j无直接联系*/ tP=*(V+j)*(Pji); tP=*(V+i)*Yji.G-tP; tP*=*(V+i); tQ=-*(V+j)*(Qji); tQ-=*(V+i)*(Yji.B-Yji.B0/Yji.k); tQ*=*(V+i); else if(Yij.B=0) conti

21、nue; /*若Bij=0，则节点i、j无直接联系*/ 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); fclose(out); system(cmd /c start out.txt); return(0);/* 主函数 */* 主函数 */* 子函数：读节点数据 */void in_node() struct Nodet

22、ype *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=(struct Nodetype *)malloc(sizeof(struct Nodetype)*No

23、de);/*节点功率*/ j=1; while(k!=0) /*若k=0，表明节点数据读取完毕*/ switch(k) case 1:k=NQ;NQ+;break; /*NQ统计PQ节点个数*/ case 2:k=NP-j;j+;break; /*从NP-1个空间倒着保存PV节点*/ case 3:k=NP;break; /*平衡节点*/ default:ErrorMsg(3); Noi=No+k; fscanf(in,%d %f %f %f %f,&i,&(*Noi).P,&(*Noi).Q,V+k,Dlta+k); i-; /*节点编号减一，以和程序表达方式兼容*/ *(Num+i)=k;

24、 /*第i个Num元素中存放i节点在No中的下标*/ fscanf(in,%d,&k); /*读取节点类型*/ if(NQ+j!=Node) ErrorMsg(4); /*检验节点数据个数是否够Node个*/fprintf(out,【 节点参数表 】n);PrtNode();fprintf(out, 总节点：%dn PQ节点：%dn PV节点：%dn,Node,NQ,NP-NQ);/* 子函数：读线路数据，并形成节点导纳矩阵Y */void in_line() struct Linetype *Yi; float R,X,k,B; m=sizeof(struct Linetype); Y=(s

25、truct Linetype *)malloc(m*Node); /*先开Node行,每一个节点一行*/ for(i=0;iNode;i+) /*再在第i行上面开辟Node-i列*/ /*即以上三角存储 Y 矩阵*/ *(Y+i)=(struct Linetype *)malloc(m*(Node-i); Yi=*(Y+i)-i; for(j=i;jj) tmp=i;i=j;j=tmp;k=1/k;B*=k; Yi=*(Y+i)-i; /* 以Yi代替*(Y+i)-i，简化表达式并避免重复计算 */ Yij.B0=B; /*保存ij0、ji0对地充电电容到Bij0*/ Yij.k=k; /*且

26、有B0ji=B0ij/k*/ tmp=R*R+X*X; R/=tmp; X/=tmp; Yij.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=0;kNode;k+) i=k; i=*(Num+i); /*查取第i节点在程序中存储序号*/ for(j=0;jk;j+) fprintf(out,ttt); for(m=k;mNode;m+) j=*(Num+m); /*查取第j节点在程序中存储序号*/ if(ij) fprintf(out,(%10.6f,%10.6f) ,Y_(i,j).G,Y_(i,j).B); else fprintf(out,(%10.6f,%10.6f) ,Y_(j,i).G,Y_(j,i).B); fprintf(out,n); /* 子函数：生成BP、BQ矩阵 */void B_Form() f