电力系统潮流计算程序实现C语言版（动态节点+稀疏技术）

概述

1. 电力系统潮流计算程序采用的是牛顿-拉夫逊法(直角坐标)

2. 潮流计算程序的系统节点数是可以动态变化的，根据节点大小分配储存空间

3. 为了减少储存空间和加快计算速度，根据电力系统节点导纳矩阵稀疏的特点，计算程序运用了稀疏储存技术，详细讲解请见参考目录。

4. 关于本程序的计算原理和流程图均可参考后面参考目录中的图书

源数据格式说明

源数据有功功率、无功功率、电压、阻抗、感抗、对地充电容纳均以标幺值表示。
数据文件必须命名为data.txt且与潮流计算程序放置于同一个文件中。数据文件data.txt包含两类参数：节点参数和线路参数。节点数据块与线路数据块之间用数字0作为间隔，即在节点数据块结束后，另起一行输入0，然后再在后面按格式要求录入线路参数。

1. 节点参数
节点参数包括：系统节点数Node、节点功率(有功P、无功Q)，节点电压(有效值V、相角Delta) ，参数组织格式：

1. 节点数Node 节点数Node写在参数文件的开头，如：4表明为四节点系统。
2. 功率和电压P/Q/V/Delta
   首先给出节点参数示例： 2 3 0.5 0 1.10 0 第一列数字2表明该节点的类型为2-PV节点；第二列数字3表明该行数据为节点3的参数；后面三列依次为P、Q、V的给定值，给定值为0，表明该项参数未知；第六列为相角δ，非平衡节点的δ即为PQ迭代的初始相角值，平衡节点的即为给定的相角值。 节点功率为各支路输入功率之和，且规定功率流入节点为正，流出为负。故而负荷功率为负值，发电机功率为正值。
3. 节点类型的判断 根据节点的给定参数可以将节点分为三种类型：
   1）PQ节点：给定P、Q初始值的节点，用数字代码1表示；
   2）PV节点：给定P、V初始值的节点，用数字代码2表示；
   3）平衡节点：没有给定P、Q初始值，仅给定V的初始值，用数字代码3表示；
   注意：数据文件中，节点顺序为PQ、PV、平衡节点

2.线路参数
线路参数包括：线路阻抗(R,X)、变压器变比k和对地充电(接地支路)容纳b。每一条线路包括节点号在内共有6个参数，6个参数缺一不可。

1. 线路阻抗Z
   首先给出线路参数示例： 1 4 0.12 0.5 1 0.01920 前两列代表线路两端的节点编号，亦即线路编号14(即41)；第三四列代表线路阻抗Z=0.12+j0.5；第五列为变压器变比k；最后一列为线路的对地充电电容的一半，即B/2。 k不为0。当k=1则表明该支路为普通支路；否则该支路为变压器支路。普通支路没有任何特殊要求，但对于变压器支路，有以下注意事项。
2. 变压器支路
   首先给出变压器线路示例： 3 1 0 0.3 0.909091 0
   1）变比k为转换为变压器支路的标准等值电路(如下图)后的变比。
   
   2）变压器线路的编号31有特定含义：3对应节点p，1对应节点q。即变压器支路的编号对节点顺序有要求，p节点编号在前，q节点在后。
   3）q节点为全系统参考电压侧。
   程序智能化读取线路数据，用户无需刻意对线路参数进行归类、排序，可随意输入，只要每行数据都各自符合格式要求即可。

输出数据说明

1. 节点参数表
节点参数表实际上只是对源数据的重现，用于对源数据的整理以及检验程序是否成功读取原始数据。 其中，节点类型列中，BS代表平衡节点、PQ代表PQ节点、PV代表PV节点。 参数表后紧随的是对各类节点数目的统计结果。 该表属于对原始数据的基本统计处理。

2. 节点导纳矩阵
节点导纳矩阵以上三角方式按照节点编号顺序输出，括号内部逗号前为导纳Gij、逗号后为容纳Bij。下三角部分可由上三角部分对称得到。

3. 潮流计算结果(节点)
该部分计算结果包含PQ迭代次数k和节点参数计算结果。迭代次数k在65535次以内均可正常显示。一般PQ迭代不会达到这个数值，否则可以认为PQ迭代是发散的。
节点参数计算结果表，展示了PQ迭代收敛时的计算结果，包括非平衡节点的电压相角δ、PQ节点的电压V、PV节点的无功功率Q和平衡节点的功率。

4. 潮流计算结果(线路)（该部分暂未编写，请等待后续完善）
线路功率计算结果表，展示了PQ迭代收敛时的线路功率计算结果S，包括线路有功功率P和无功功率Q。

参考文献

1. 陈亚明.电力系统计算程序.北京：水利电力出版社，1995
2. 张伯明，陈寿孙，严正.高等电力网络分析（第2版）.北京：清华大学出版社，2007

源代码

/*   FUNCTION : POWER FLOW             */ /*   WRITTEN BY : ZAHGN CHUN           */ /*   LAST EDITED : 2014-12-10          */ #include <stdio.h> #include <math.h> #include<stdlib.h> #include<string.h> #include<conio.h> #include<iostream>#include<vector>using namespace std;/***  子函数声明  ***/ void in_data();   /* 读取数据 */ void Y_Form();    /* 生成导纳矩阵 */void chuzhi();    /* 给节点电压赋初值 */   void B_Form();    /* 生成雅克比矩阵和求解修正方程*/  void xinzhi();    /*计算新值*/ void PrtNode();   /*输出潮流计算结果*/ void PrtY();      /*打印Y矩阵*/ /***  全局变量声明  ***/ int Node;              /*网络节点数*/ int Num;               /*网络支路数*/ int i,j,k,k1,k2,k5,l,m,ig,i1,i2,i3,i5,t,s;          /*通用下标值*/ int J1,J0,J3;float tmp;            /*临时数据暂存*/ FILE *in,*out;        /*输入、输出文件指针*/ int  J5;                //雅克比矩阵第i行非零元素的技术单元int  k0=1;             //导纳矩阵非零互导纳元素的指针int  FL=0;             //导纳矩阵非零互导纳元素的个数float   a1,b1,r1,x1,a,r,x,b;       //寄存单元float   gij,bij,gi,bi,gj,bj;float A3;            //第i节点电流的实部float B3;            //第i节点电流的虚部float P2,Q2,P3,Q3,C3,D3,C5,D5,A5,B5,C6;float  PG,QG,PL,QL;     //全网功率float  E,F,PI;int N0;               //平衡节点的节点号float U0;           //存平衡节点的给定电压值int IQ;               //发电机台数int IP;               //负荷个数int N1;               //PV节点个数float UP=1;;             //PQ节点的电压初值int   it;              //迭代次数float  am=0;             //节点功率误差的最大绝对值int    i0;             //节点功率误差绝对值最大的节点号//声明变长数组vector<int> IZA(Num);          //存支路状态数vector<int> IZ1(Num);          //存支路一端的节点号vector<int> IZ2(Num);          //存支路另一端的节点号vector<float> Z1(Num);          //存支路正序电阻vector<float> Z2(Num);          //存支路正序电抗vector<float> Z3(Num);          //存支路正序电纳或非标准变比vector<int> IWGA(Node);          //存发电机状态数vector<int> IWG(Node);          //存发电机节点号vector<float> W1(Node);          //存发电机有功功率vector<float> W2(Node);          //存发电机无功功率vector<int> ILP(Node);          //存负荷状态数vector<int> ILD(Node);          //存负荷节点号vector<float> WL1(Node);          //存负荷有功功率vector<float> WL2(Node);          //存负荷无功功率vector<int> IPV(5);             //存PV节点的节点号vector<float> PV(5);            //存PV节点的给定电压值vector<float> YZ1(Node*Node);     //不规则互导纳元素的实部vector<float> YZ2(Node*Node);     //不规则互导纳元素的虚部vector<int> IY1(Node*Node);            //不规则互导纳元素的行号vector<int> IY2(Node*Node);            //不规则互导纳元素的列号vector<float> D1(Node);           //导纳矩阵对角元素的实部    vector<float> D2(Node);           //导纳矩阵对角元素的虚部vector<float> Y1(Node*Node);     //导纳矩阵非对角非零元素的实部vector<float> Y2(Node*Node);     //导纳矩阵非对角非零元素的虚部vector<int> IY(Node*Node);       //导纳矩阵非对角非零元素的列号vector<int> IN(Node);            //导纳矩阵每行非对角非零元素的个数vector<float> U1(Node);           //节点电压的实部vector<float> U2(Node);           //节点电压的虚部vector<float> AK1(Node*Node);           //存雅克比矩阵一行的非零元素vector<float> AK2(Node*Node);vector<float> AK3(Node*Node);vector<float> AK4(Node*Node);vector<int>   JK(Node*Node);            //存雅克比矩阵一行的非零元素的列号vector<float> AJ1(Node*Node);        //存雅克比矩阵消去规格化后上三角非对角非零元素vector<float> AJ2(Node*Node);vector<float> AJ3(Node*Node);vector<float> AJ4(Node*Node);vector<int> IJ(Node*Node);            //存雅克比矩阵消去规格化后上三角非对角非零元素的列号vector<int> JF(Node);             //存雅克比矩阵消去规格化后上三角每行非对角非零元素的个数vector<float> CK1(Node);          //先存节点无功功率误差ΔQi（或ΔUi的平方），后存节点电压修正量的实部vector<float> CK2(Node);          //先存节点有功功率误差ΔPi，后存节点电压修正量的虚部vector<float> PD(Node);           //节点负荷的有功功率vector<float> QD(Node);           //节点负荷的无功功率vector<float> PF(Node);           //节点给定的发电有功功率vector<float> QF(Node);           //节点给定的发电无功功率vector<float> GP(Node);           //节点实际的发电有功功率vector<float> GQ(Node);           //节点实际的发电无功功率vector<float> Q(Node);           //节点的无功功率误差值vector<int> IVI(Node);             //PV节点的标致/******************  主函数  ******************/ /**↓****↓****↓**  主函数  **↓****↓****↓**/ int main(void) {   if((in=fopen("Data.txt","r"))==NULL) printf("wrong\n");   if((out=fopen("out.txt","w"))==NULL) printf("wrong\n");   in_data();   /* 读取数据 */    for(i=0;i<Node*Node;i++)  {      Y1[i]=0; Y2[i]=0; IY[i]=0; YZ1[i]=0; YZ2[i]=0;IY1[i]=0;      AK1[i]=AK2[i]=AK3[i]=AK4[i]=0;      AJ1[i]=AJ2[i]=AJ3[i]=AJ4[i]=0;      JK[i]=JF[i]=IJ[i]=0;  }  for(i=0;i<Node;i++)  {   D1[i]=D2[i]=0;   }  Y_Form();    PrtY();  chuzhi();  for(it=1;it<20;it++)  {       B_Form();        fprintf(out,"第%d迭代节点功率误差绝对值最大的节点号:%d,误差的最大绝对值:%f\n",it,i0,am);       for(i=0;i<Node-1;i++)       fprintf(out,"e[%d]=%f f[%d]=%f\n",i,CK1[i],i,CK2[i]);     printf("IT=%d\n",it);       if(am<0.0001)          break;       xinzhi();  }  PrtNode();      /*输出潮流计算结果*/   fclose(out);   system("cmd /c start out.txt");   return(0); } /**↑****↑****↑**  主函数  **↑****↑****↑**/ /**************** 计算新值  ****************/ void xinzhi(){        for(i=Node-2;i>=0;i--)                 //解修正方程回代运算             //zhuyi!!!!!!!!!!!!!!!!    {                       for(ig=0;ig<JF[i];ig++)        {            k--;            m=IJ[k];            CK1[i]=CK1[i]-AJ1[k]*CK1[m]-AJ2[k]*CK2[m];            CK2[i]=CK2[i]-AJ3[k]*CK1[m]-AJ4[k]*CK2[m];        }    }    for(i=0;i<Node;i++)               //计算节点新的电压值    {        U1[i]=U1[i]+CK1[i];        U2[i]=U2[i]+CK2[i];    }    for(i=0;i<Node-1;i++)       printf("e[%d]=%f f[%d]=%f\n",i,U1[i],i,U2[i]);}/******** 生成雅克比矩阵和不平衡量  *********/void  B_Form(){    am=0;     k0=0;    for(i=0;i<Node;i++)    {                           a=D1[i]*U1[i];        b=D2[i]*U1[i];        r=D1[i]*U2[i];        x=D2[i]*U2[i];        A3=a-x;        B3=r+b;        J5=1;                            for(ig=0;ig<IN[i];ig++)         {            j=IY[k0];            A3=A3+Y1[k0]*U1[j]-Y2[k0]*U2[j];            B3=B3+Y1[k0]*U2[j]+Y2[k0]*U1[j];            if(i!=N0&&j!=N0)                                {                JK[J5]=IY[k0];                                                        AK1[J5]=Y1[k0]*U2[i]-Y2[k0]*U1[i];                AK3[J5]=Y1[k0]*U1[i]+Y2[k0]*U2[i];                AK2[J5]=-AK3[J5];                AK4[J5]=AK1[J5];                       J5++;            }            k0++;        }        if(i==N0)        {            GQ[i]=A3*U2[i]-B3*U1[i];            GP[i]=A3*U1[i]+B3*U2[i];            CK1[i]=0;            CK2[i]=0;            JF[i]=0;        }        else        {            P2=PD[i]+(A3*U1[i]+B3*U2[i]);            Q2=QD[i]+(A3*U2[i]-B3*U1[i]);            GP[i]=P2;            GQ[i]=Q2;            P2=PF[i]-P2;            Q2=QF[i]-Q2;            P3=P2;            Q3=Q2;            if(IVI[i]==1)                if(Q[i]>=0||it-1<5)                                                  Q3=0;            AK3[0]=A3+a+x;            AK4[0]=B3-b+r;            JK[0]=i;                                CK2[i]=P2;            if(IVI[i]==1)            {                Q[i]=Q2;                if(it>=5)                {                    if(Q2>=0)                    {                        for(j=1;j<J5;j++)                                             AK1[j]=AK2[j]=0;                          AK1[0]=2*U1[i];                        AK2[0]=2*U2[i];                        Q2=0;                        CK1[i]=Q2;                    }                    else                    {                        AK1[0]=r-b-B3;                        AK2[0]=A3-a-x;                        CK1[i]=Q2;                    }                }                else                {                        for(j=1;j<J5;j++)                                              AK1[j]=AK2[j]=0;                          AK1[0]=2*U1[i];                        AK2[0]=2*U2[i];                        Q2=0;                        CK1[i]=Q2;                }            }            else            {                AK1[0]=r-b-B3;                AK2[0]=A3-a-x;                CK1[i]=Q2;            }            C5=fabs(P3);            D5=fabs(Q3);            if(C5>am)            {   i0=i;am=C5;   }            if(D5>am)            {   i0=i;am=D5;   }            k=0;                               //消去运算            m=0;            for(i1=0;i1<i;i1++)            {                      for(i3=1;i3<J5;i3++)                    if(JK[i3]==i1)  break;                if(i3<J5)                {      //  printf("xiaoq\n");                      for(ig=0;ig<JF[i1];ig++)                    {                                           for(i2=0;i2<J5;i2++)                            if(JK[i2]==IJ[k])   break;                        if(i2<J5)                        {                                             AK1[i2]=AK1[i2]-AK1[i3]*AJ1[k]-AK2[i3]*AJ3[k];                                      AK2[i2]=AK2[i2]-AK1[i3]*AJ2[k]-AK2[i3]*AJ4[k];                            AK3[i2]=AK3[i2]-AK3[i3]*AJ1[k]-AK4[i3]*AJ3[k];                            AK4[i2]=AK4[i2]-AK3[i3]*AJ2[k]-AK4[i3]*AJ4[k];                         }                        else                        {                               AK1[J5]=-AK1[i3]*AJ1[k]-AK2[i3]*AJ3[k];                            AK2[J5]=-AK1[i3]*AJ2[k]-AK2[i3]*AJ4[k];                            AK3[J5]=-AK3[i3]*AJ1[k]-AK4[i3]*AJ3[k];                            AK4[J5]=-AK3[i3]*AJ2[k]-AK4[i3]*AJ4[k];                                            JK[J5]=IJ[k];                            J5++;                        }                        k++;                            }                    CK1[i]=CK1[i]-AK1[i3]*CK1[i1]-AK2[i3]*CK2[i1];                    CK2[i]=CK2[i]-AK3[i3]*CK1[i1]-AK4[i3]*CK2[i1];                }                else                {                    k=k+JF[i1];                    }            }            k5=k;                     //规格化            a=AK1[0]*AK4[0]-AK2[0]*AK3[0];                    if(a==0)            {                CK1[i]=0;                CK2[i]=0;                JF[i]=0;            }            else            {                A3=AK4[0]/a;                B3=-AK2[0]/a;                C3=-AK3[0]/a;                D3=AK1[0]/a;                for(j=1;j<J5;j++)                    if(JK[j]>i)                    {                                               //  printf("guife\n");                         IJ[k]=JK[j];                                   AJ1[k]=A3*AK1[j]+B3*AK3[j];                        AJ2[k]=A3*AK2[j]+B3*AK4[j];                        AJ3[k]=C3*AK1[j]+D3*AK3[j];                        AJ4[k]=C3*AK2[j]+D3*AK4[j];                        k++;                    }                A5=CK1[i];                B5=CK2[i];                CK1[i]=A3*A5+B3*B5;                CK2[i]=C3*A5+D3*B5;                            JF[i]=k-k5;            }          }    } }/****************  子函数：读数据  ****************/ void in_data() {     fscanf(in,"%d %d %d %d %d",&Node,&Num,&IQ,&IP,&N1);       /*读取节点数Node、支路数、发电机个数、负荷个数、PV节点个数*/     for(i=0;i<Num;i++)    {     Z1[i]=Z2[i]=Z3[i]=0;  }    for(i=0;i<IQ;i++)    {     IWGA[i]=IWG[i]=W1[i]=W2[i]=0; }   for(i=0;i<Num;i++)      fscanf(in,"%d %d %d %f %f %f",&IZA[i],&IZ1[i],&IZ2[i],&Z1[i],&Z2[i],&Z3[i]);    for(i=0;i<IQ;i++)       fscanf(in,"%d %d %f %f",&IWGA[i],&IWG[i],&W1[i],&W2[i]);    for(i=0;i<IP;i++)   {       fscanf(in,"%d %d %f %f",&ILP[i],&ILD[i],&WL1[i],&WL2[i]);    }   for(i=0;i<N1;i++)     fscanf(in,"%d %f",&IPV[i],&PV[i]);     fscanf(in,"%d %f",&N0,&U0); } /****************给节点电压赋初值****************/ void chuzhi()    {    for(i=0;i<Node;i++)    {        PD[i]=0;        QD[i]=0;        PF[i]=0;        QF[i]=0;        IVI[i]=0;    }    for(i=0;i<IP;i++)    {        if(ILP[i]!=0)        {            j=ILD[i];            PD[j]=WL1[i];            QD[j]=WL2[i];        }    }    for(i=0;i<IQ;i++)    {        if(IWGA[i]!=0)        {            j=IWG[i];            PF[j]=W1[i];            QF[j]=W2[i];        }    }    for(i=0;i<N1;i++)    {        j=IPV[i];        IVI[j]=1;    }    for(i=0;i<Node;i++)    {        U1[i]=UP;          U2[i]=0;    }    for(i=0;i<N1;i++)    {        j=IPV[i];              U1[j]=PV[i];    }    U1[N0]=U0;}/************形成节点导纳矩阵Y  ************/ void Y_Form(){    for(i=0;i<Node;i++)    {        D1[i]=0;        D2[i]=0;    }    l=0;                                for(k=0;k<Num;k++)    {        ig=IZA[k];        i=IZ1[k];                    j=IZ2[k];          r=Z1[k];            x=Z2[k];        b=Z3[k];        a=r*r+x*x;          if(a!=0)        {            gij=r/a;bij=-x/a;                      if(ig==1)            {                gi=gij;                gj=gij;                bi=bij+b;                            bj=bi;                D1[i]=D1[i]+gi;                          D2[i]=D2[i]+bi;                          D1[j]=D1[j]+gj;                          D2[j]=D2[j]+bj;                          YZ1[l]=-gij;                                              YZ2[l]=-bij;                                      IY1[l]=i;                   IY2[l]=j;                    l++;                YZ1[l]=-gij;                                     YZ2[l]=-bij;                                     IY1[l]=j;                IY2[l]=i;                    l++;            }            else            {                if(ig==2||ig==3)                {                    gj=gij;                    bj=bij;                    gij=gij/b;                    bij=bij/b;                    gi=gij/b;                    bi=bij/b;                    D1[i]=D1[i]+gi;                         D2[i]=D2[i]+bi;                        D1[j]=D1[j]+gj;                        D2[j]=D2[j]+bj;                       YZ1[l]=-gij;                                YZ2[l]=-bij;                               IY1[l]=i;                                  IY2[l]=j;                      l++;                    YZ1[l]=-gij;                               YZ2[l]=-bij;                                IY1[l]=j;                    IY2[l]=i;                       l++;                }                else                {                    D1[i]=D1[i]+gij;                    D2[i]=D2[i]+bij;                }            }        }    }    j=0;    k0=0;    m=0;     //出口标志    for(i=0;i<Node;i++)    {        J1=0;        for(k=0;k<l;k++)        {                                    if(IY1[k]==i)            {                J3=IY2[k];                for(k1=0;k1<J1;k1++)                {                    k2=k0+k1;                                 if(J3==IY[k2])                                    {   break;  m=1;  }                }                if(m==1)                {                    Y1[k2]=Y1[k2]+YZ1[k];                          Y2[k2]=Y2[k2]+YZ2[k];                }                else                {                                     Y1[j]=YZ1[k];                    Y2[j]=YZ2[k];                    IY[j]=J3;                          J1++;                    j++;                }            }        }        IN[i]=J1;        k0=k0+J1;    }    FL=k0;}/*******************打印导纳矩阵**************/ void PrtY(){     fprintf(out,"非对角线元素：\n");      for(i=0;i<FL;i++)      {        fprintf(out,"%10.6f,%10.6f,,,,,,%d\n",Y1[i],Y2[i],IY[i]);     }     fprintf(out,"\n");      fprintf(out,"对角线元素：\n");      for(i=0;i<Node;i++)      {        fprintf(out,"%10.6f,%10.6f\n",D1[i],D2[i]);     }     fprintf(out,"\n");     fprintf(out,"每行非对角线非零元素个数：\n");      for(i=0;i<Node;i++)      {         fprintf(out,"%d\n",IN[i]);     }     fprintf(out,"\n"); }/***************输出潮流计算结果****************/void PrtNode(){    GP[N0]=GP[N0]+PD[N0];              //输出节点信息    GQ[N0]=GQ[N0]+QD[N0];    PG=0;    QG=0;    PL=0;    QL=0;    PI=180/3.14159;    fprintf(out,"\n");    for(i=0;i<Node;i++)    {        E=U1[i];        F=U2[i];        a=sqrt(E*E+F*F);        b=PI*atan2(F,E);        PG=PG+GP[i];        QG=QG+GQ[i];        PL=PL+PD[i];        QL=QL+QD[i];        fprintf(out,"第%d节点：\n幅值：%f，相角：%f\n",i+1,a,b);    } // fprintf(out,"全网发电有功功率：%f，无功功率：%f，负荷有功功率：%f，负荷无功功率：%f\n",PG,QG,PL,QL);       //输出支路信息    float ploss=0,qloss=0;       //全网有功、无功功率损耗    float qb=0;                //全网输电线充电功率}