动态节点

概述

  1. 电力系统时髦统计程序行使的是牛顿(Newton)-拉夫逊法(直角坐标)

  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为转移为变压器支路的业内等值电路(如下图)后的变比。
    图片 1
    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;                //全网输电线充电功率

}