附录_潮流计算程序

附录  牛顿-拉夫逊法潮流计算程序

本附录所介绍的牛顿-拉夫逊法潮流计算程序采用极坐标形式，其中所涉及的计算公式与第四章中的式(4-42)-(4-59)完全相同，计算程序框图与图4-6基本一致。程序采用C语言编制。为了便于初学者阅读，节点导纳矩阵和雅可比矩阵都用满阵存储而未采用稀疏技巧；但为了适当照顾使用的方便性，在输入数据中对节点编号次序不作任何要求。

下面先介绍输入文件的格式和要求，然后列出程序，最后说明潮流结果的输出。建议读者先从原始数据的输入中了解和熟悉它们在程序中对应的变量、结构体数组及其成员的名称，然后对照图4-6和式(4-42)-(4-59)仔细和耐心地阅读程序，最好能在计算机上亲自实现和进行调试。在调试时，可以用下面给出的对应于[例4-3]系统的输入数据，以及程序运行中得出的中间结果，逐步与[例4-3]中所给出的中间结果进行对比，从而查找错误所在并进行改正。

一、原始数据的输入

程序通过“输入数据.txt”文件输入以下5个数据段。

1.信息(共6个)

(1)总节点数(变量num_node)

(2)线路和并联电容器总数(变量num_line)

(3)变压器支路总数(变量num_tran)

(4)发电机节点总数(变量num_gene)

(5)负荷节点总数(变量num_load)

(6)节点功率不平衡量的容许误差(变量error)

2.线路和并联电容器数据(结构体数组line)：每一线路或并联电容器包括5个数据

线路	并联电容器	成员
I侧节点编号	所接节点编号	i
J侧节点编号	同上	j
Π型等值电路电阻	电容器电阻	a
Π型等值电路电抗	电容器电抗(负数)	b
Π型等值电路一端电纳	0.0	c

I侧节点编号和J侧节点编号可以对换；线路和并联电容器之间的次序可以任意，而且允许多条线路或多个电容器并联。

3.变压器支路数据(结构体数组tran)：每一变压器支路包括5个数据

变压器支路	成员
1侧节点编号	i
2侧节点编号	j
电阻	a
电抗	b
非标准变比	c

变压器电阻、电抗和非标准变比与两侧节点编号之间的关系服从图2-29,即电阻和电抗在1侧而非标准变比在2侧；三绕组变压器需按图2-26化成3个变压器支路，其中3侧变压器支路的非标准变比为1。变压器支路之间的次序可以任意。

4.发电机节点数据(结构体数组gene)：每一发电机节点包括5个数据

发电机节点	成员
所在节点编号	i
节点种类	j
发出有功功率	a
发出无功功率	b
电压	c

对于PQ节点，节点种类为1，电压可给任意值；对于PV节点，节点种类为-1，发出无功功率可给任意值；对于平衡节点，节点种类为0，发出有功功率和发出无功功率都可给任意值。发电机节点之间的次序可以任意。

5.负荷节点数据(结构体数组load)：每一负荷节点包括3个数据

负荷节点	成员
所在节点编号	i
吸收有功功率	a
吸收无功功率	b

负荷机节点之间的次序可以任意。

对于[例4-3]中的电力系统，输入数据如下：

4   4   1   2   2   0.00001

4  3  0.260331     0.495868    0.0258

1  4  0.173554     0.330579    0.017243

2  2  0.000000   -20.000000    0.000000    

3  1  0.130165     0.247934    0.012932

1  2  0.000000     0.166667    1.128205

4   0  0.0   0.0   1.05

3  -1  0.2   0.0   1.05

2  0.5   0.3

4  0.15  0.1

二、计算程序

    /* * * * * *     牛顿-拉夫逊法潮流计算程序     * * * * * */

# include

# include

# include

# define pnt         1     // 1 - 输出中间结果

void read_data();          // 输入始数据

void admt_matrix();        // 形成导纳矩阵

void form_Jacobian();      // 形成雅可比矩阵, 计算功率误差

void solv_Eqn();           // 求解修正方程式

void node_flow();          // 输出节点潮流

void branch_flow();        // 输出支路潮流

double **newSpaceDouble2(int,int);

void deleteSpaceDouble2(double **,int);

int num_node,num_line,num_tran,num_gene,num_load,iter;

struct data *line,*tran,*gene,*load;

double **G,**B,**Jacob;

double *Um,*Ua,*P,*Q;

double error_max;

FILE *fin,*fou,*chk;

struct data

{

    int i;

    int j;

    double a;

    double b;

    double c;

};

void main()

{

    int i,j,conv;

    double a,error;

    fin=fopen("输入数据.txt

    if(fin==NULL)

    {

        printf(" 注意! 没有“输入数据.txt”文件\n");

        exit(0);

    }

    fou=fopen("潮流输出.txt

    if(pnt==1) chk=fopen("中间结果.txt

// 输入原始数据和形成节点导纳矩阵

    read_data();

    G=newSpaceDouble2(num_node,num_node);

    B=newSpaceDouble2(num_node,num_node);

for(i=1;i<=num_node;i++)

        for(j=1;j<=num_node;j++)

            G[i][j]=B[i][j]=0.0;

    admt_matrix();

// 给定电压有效值和相位初值

    Um=new double[num_node+1];

    Ua=new double[num_node+1];

for(i=1;i<=num_node;i++)

    {

        Um[i]=1.0;

        Ua[i]=0.0;

    }

for(i=1;i<=num_gene;i++)

        if(gene[i].j<=0) Um[gene[i].i]=gene[i].c;

    iter=0;

// 形成雅可比矩阵计算功率误差

    Jacob=newSpaceDouble2(2*num_node,2*num_node+1);

    P=new double[num_node+1];

    Q=new double[num_node+1];

R2:    form_Jacobian();

// 收敛判断

    error=0.0;

for(i=1;i<=2*num_node;i++)

        if(fabs(Jacob[i][2*num_node+1])>error)

            error=fabs(Jacob[i][2*num_node+1]);

    fprintf(fou,"\\n      迭代次数: %2d        最大功率误差: %11.6f

        iter+1,error);

    if(error    {

        conv=1;

        goto R1;

    }

    if((iter>10) || (error>1.0e4)) // 潮流计算不收敛

    {

        fprintf(fou,"\\n\\n      潮流不收敛");

        goto nd;

    }

// 求解修正方程式并修正电压

    solv_Eqn();

for(i=1;i<=num_node;i++)

    {

        a=Jacob[i][2*num_node+1];

        Ua[i]=Ua[i]+(-1*a);

        a=Jacob[num_node+i][2*num_node+1];

        Um[i]=Um[i]-(Um[i]*a);

    }

if((pnt==1) && (iter<2))

    {

        fprintf(chk,"\\n\\n  电压相位和有效值新值, 迭代 %d\\n\\n",iter+1);

        for(i=1;i<=num_node;i++)

            fprintf(chk," %3d  %10.5f  %8.5f \\n",i,Ua[i],Um[i]);

    }

    iter=iter+1;

    goto R2;

// 输出潮流结果

R1:    node_flow();

    branch_flow();

nd:    fclose(fin);

    fclose(fou);

    if(pnt==1) fclose(chk);

    free(line); free(tran); free(gene); free(load);

    deleteSpaceDouble2(G,num_node);

    deleteSpaceDouble2(B,num_node);

    deleteSpaceDouble2(Jacob,2*num_node);

    free(Um); free(Ua); free(P); free(Q);

}

void read_data() // 输入始数据

{

    int i;

    fscanf(fin,"%d %d %d %d %d %lf",&num_node,&num_line,

        &num_tran,&num_gene,&num_load,&error_max);

    line=(struct data *)calloc(num_line+1,sizeof(struct data));

    tran=(struct data *)calloc(num_tran+1,sizeof(struct data));

    gene=(struct data *)calloc(num_gene+1,sizeof(struct data));

    load=(struct data *)calloc(num_load+1,sizeof(struct data));

for(i=1;i<=num_line;i++) fscanf(fin,"%d %d %lf %lf %lf

       &line[i].i,&line[i].j,&line[i].a,&line[i].b,&line[i].c);

for(i=1;i<=num_tran;i++) fscanf(fin,"%d %d %lf %lf %lf

        &tran[i].i,&tran[i].j,&tran[i].a,&tran[i].b,&tran[i].c);

for(i=1;i<=num_gene;i++) fscanf(fin,"%d %d %lf %lf %lf

        &gene[i].i,&gene[i].j,&gene[i].a,&gene[i].b,&gene[i].c);

for(i=1;i<=num_load;i++) fscanf(fin,"%d%lf %lf

        &load[i].i,&load[i].a,&load[i].b);

}

void admt_matrix() // 形成导纳矩阵

{

    int i,j;

    double r,x,b,kt;

    struct data *p,*end;

    end=line+num_line;

    for(p=line+1;p<=end;p++) // 线路

    {

        i=p->i;

        j=p->j;

        r=p->a;

        x=p->b;

        b=r*r+x*x;

        r=r/b;

        x=-x/b;

        if(i==j)

        {

            G[i][i]+=r;

            B[i][i]+=x;

            continue;

        }

        b=p->c;

        G[i][j]=G[i][j]-r;

        B[i][j]=B[i][j]-x;

        G[j][i]=G[i][j];

        B[j][i]=B[i][j];

        G[i][i]=G[i][i]+r;

        B[i][i]=B[i][i]+x+b;

        G[j][j]=G[j][j]+r;

        B[j][j]=B[j][j]+x+b;

    }

    end=tran+num_tran; // 变压器

    for(p=tran+1;p<=end;p++)

    {

        i=p->i;

        j=p->j;

        r=p->a;

        x=p->b;

        b=r*r+x*x;

        r=r/b;

        x=-x/b;

        kt=p->c;

        G[i][i]+=r;

        B[i][i]+=x;

        G[i][j]=G[i][j]-r/kt;

        B[i][j]=B[i][j]-x/kt;

        G[j][i]=G[i][j];

        B[j][i]=B[i][j];

        r=r/kt/kt;

        x=x/kt/kt;

        G[j][j]+=r;

        B[j][j]+=x;

    }

    if(pnt==1)

    {

        fprintf(chk,"\\n\\n            导纳矩阵中的非零元素\n\\n");

        for(i=1;i<=num_node;i++) for(j=i;j<=num_node;j++) 

            if((G[i][j]!=0.0) || (B[i][j]!=0.0))

                fprintf(chk," %3d %4d %16.5f %16.5f\\n",i,j,G[i][j],B[i][j]); 

    }

}

void form_Jacobian() // 形成雅可比矩阵, 计算功率误差

{

    int i,j,nu,ii,k,n2,kk;

    double vi,di,dij,vj,dj,p,q,b,g,gp,gq,lp,lq;

    double Hij,Lij,Nij,Mij,Hii,Nii,Mii,Lii,dp,dq;

    nu=2*num_node+1;

    n2=2*num_node;

    for(i=1;i<=n2;i++)

        for(j=1;j<=nu;j++)

            Jacob[i][j]=0.0;

for(i=1;i<=num_node;i++)

    {

        vi=Um[i];

        di=Ua[i];

        dp=0.0;

        dq=0.0;

        for(j=1;j<=num_node;j++) // 非对角元素

        {

            if(j==i) continue;

            g=G[i][j];

            b=B[i][j];

            vj=Um[j];

            dj=Ua[j];

            dij=di-dj;

            Hij=-Um[i]*Um[j]*(g*sin(dij)-b*cos(dij));

            Lij=Hij;

            Jacob[i][j]=Hij;

            Jacob[i+num_node][j+num_node]=Lij;

            Nij=-Um[i]*Um[j]*(g*cos(dij)+b*sin(dij));

            Mij=-Nij;

            Jacob[i][j+num_node]=Nij;

            Jacob[i+num_node][j]=Mij;

            p=Um[j]*(g*cos(dij)+b*sin(dij));

            q=Um[j]*(g*sin(dij)-b*cos(dij));

            dp=dp+p;

            dq=dq+q;

        }

        g=G[i][i];

        b=B[i][i];

        Hii=vi*dq;

        Nii=-vi*dp-2*vi*vi*g;

        Mii=-vi*dp;

        Lii=-vi*dq+2*vi*vi*b;

        Jacob[i][i]=Hii;

        Jacob[i][i+num_node]=Nii;

        Jacob[i+num_node][i]=Mii;

        Jacob[i+num_node][i+num_node]=Lii;

        Jacob[i][nu]=-vi*(dp+vi*g);

        Jacob[i+num_node][nu]=-vi*(dq-vi*b);

        P[i]=vi*(dp+vi*g); // 节点注入有功

        Q[i]=vi*(dq-vi*b); // 节点注入无功

    }     

for(i=1;i<=num_load;i++)

    {

        kk=load[i].i;

        lp=load[i].a;

        lq=load[i].b;

        Jacob[kk][nu]=-lp+Jacob[kk][nu];

        Jacob[kk+num_node][nu]=-lq+Jacob[kk+num_node][nu];

    }

for(i=1;i<=num_gene;i++)

    {

        kk=gene[i].i;

        gp=gene[i].a;

        gq=gene[i].b;

        Jacob[kk][nu]=gp+Jacob[kk][nu];

        Jacob[kk+num_node][nu]=gq+Jacob[kk+num_node][nu];

    }    

    for(k=1;k<=num_gene;k++) // 去掉pv及平衡节点

    {

        ii=gene[k].i;

        kk=gene[k].j;

        if(kk==0) // 平衡节点

        {

            for(j=1;j<=n2;j++)

            {

                Jacob[ii][j]=0.0;

                Jacob[num_node+ii][j]=0.0;

                Jacob[j][ii]=0.0;

                Jacob[j][num_node+ii]=0.0;

            }

            Jacob[ii][ii]=1.0;

            Jacob[num_node+ii][num_node+ii]=1.0;

            Jacob[ii][nu]=0.0;

            Jacob[num_node+ii][nu]=0.0;

        }

        if(kk<0) // pv节点

        {

            for(j=1;j<=n2;j++)

            {

                Jacob[num_node+ii][j]=0.0;

                Jacob[j][num_node+ii]=0.0;

            }

            Jacob[num_node+ii][num_node+ii]=1.0;

            Jacob[num_node+ii][nu]=0.0;

        }

    }

if((pnt==1) && (iter<2))

    {

        fprintf(chk,"\\n\\n 雅可比矩阵中的非零元素, 迭代 %d\\n\\n",iter+1);

        for(i=1;i<=2*num_node;i++)               

            for(j=1;j<=2*num_node+1;j++) 

                if((Jacob[i][j]!=0.0) && (Jacob[i][j]!=1.0))

                    fprintf(chk," %3d  %3d  %15.5f \\n",i,j,Jacob[i][j]); 

    }

}

void solv_Eqn() // 求解修正方程式

{

    int i,j,n2,nu,i1,k;

    double d,e;

    n2=2*num_node;

    nu=n2+1;

    for(i=1;i<=n2;i++) // 消去

    {

        i1=i+1;

        d=1.0/Jacob[i][i];

        for(j=i1;j<=nu;j++)

        {

            e=Jacob[i][j];

            if(e==0.0) continue;

            Jacob[i][j]=e*d;

        }

        if(i==n2) continue;

        for(j=i1;j<=n2;j++)

        {

            e=Jacob[j][i];

            if(e==0.0) continue;

            for(k=i1;k<=nu;k++)

                Jacob[j][k]=Jacob[j][k]-Jacob[i][k]*e;

        }

    }

    for(k=2;k<=n2;k++) // 回代

    {

        i=n2-k+1;

        i1=i+1;

        for(j=i1;j<=n2;j++)

            Jacob[i][nu]=Jacob[i][nu]-Jacob[i][j]*Jacob[j][nu];

    }

if((pnt==1) && (iter<2))

    {

        fprintf(chk,"\\n\\n 电压相位和有效值修正量, 迭代 %d\\n\\n",iter+1);

        for(i=1;i<=num_node;i++)               

        fprintf(chk," %3d  %10.5f  %8.5f \\n

            i,-Jacob[i][2*num_node+1],-Jacob[i+num_node][2*num_node+1]);

    }

}

void node_flow() // 输出节点潮流

{

    int i,j,k,ii,kk;

    double b1,b2,c1,c2;

    fprintf(fou,"\\n\\n\\n                     *-*-*- 潮 流 计 算 结 果 *-*-*-");

    fprintf(fou,"\\n\\n                        ——— 节 点 潮 流 ———");

    fprintf(fou,"\\n\\n    no.i      Um         Ua         PG         QG");

    fprintf(fou,"         PL         QL\\n\\n");

for(i=1;i<=num_node;i++)

    {

        b1=b2=c1=c2=0.0;

        for(j=1;j<=num_gene;j++)

        {

            ii=gene[j].i;

            kk=gene[j].j;

            if((i==ii) && (kk==0))  // 平衡节点

            {

                b1=P[ii];

                b2=Q[ii];

                for(k=1;k<=num_load;k++)

                {

                    ii=load[k].i;

                    if(i==ii)

                    {

                        c1=load[k].a;

                        c2=load[k].b;

                        b1=b1+c1;

                        b2=b2+c2;

                    }

                }

                break;

            }

            if((i==ii) && (kk==-1)) // pv节点

            {

                b1=gene[j].a;

                b2=Q[ii];

                for(k=1;k<=num_load;k++)

                {

                    ii=load[k].i;

                    if(i==ii)

                    {

                        c1=load[k].a;

                        c2=load[k].b;

                        b2=b2+c2;

                    }

                }

                break;

            }

        }

        for(j=1;j<=num_load;j++)

        {

            ii=load[j].i;

            if(i==ii)

            {

                c1=load[j].a;

                c2=load[j].b;

                break;

            }

        }

        fprintf(fou," %6d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f\\n

            i,Um[i],Ua[i]*180.0/3.141592654,b1,b2,c1,c2);

    }

}

void brabch_flow() // 输出支路潮流

{

    int i,j;

    double r,x,t,b,dij,cd,sd,ri,rj,xi,xj;

    double vi,vj,vij,pij,qij,pji,qji,dpb,dqb,ph,qh;

    struct data *p,*end;

    fprintf(fou,"\\n\\n                        ——— 支 路 潮 流 ———");

    fprintf(fou,"\\n\\n   i    j       Pij        Qij        Pji");

    fprintf(fou,"        Qji        dP         dQ\\n\\n");

    ph=0.0;

    qh=0.0;

for(end=line+num_line,p=line+1;p<=end;p++)

    {

        i=p->i;

        j=p->j;

        r=p->a;

        x=p->b;

        b=r*r+x*x;

        if(i==j)

        {

            vi=Um[i];

            b=vi*vi/b;

            pij=r*b;

            qij=x*b;

            pji=0.0;

            qji=0.0;

            dpb=pij;

            ph=ph+dpb;

            dqb=qij;

            qh=qh+dqb;

        }

        else

        {

            r=r/b;

            x=-x/b;

            b=p->c;

            dij=Ua[i]-Ua[j];

            vi=Um[i];

            vj=Um[j];

            vij=vi*vj;

            vi=vi*vi;

            vj=vj*vj;

            cd=vij*cos(dij);

            sd=vij*sin(dij);

            pij=vi*r-r*cd-x*sd;

            pji=vj*r-r*cd+x*sd;

            dpb=pij+pji;

            ph=ph+dpb;

            qij=-vi*(b+x)+x*cd-r*sd;

            qji=-vj*(b+x)+x*cd+r*sd;

            dqb=qij+qji;

            qh=qh+dqb;

        }

        fprintf(fou," %3d  %3d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f\\n

            i,j,pij,qij,pji,qji,dpb,dqb);

    }

for(end=tran+num_tran,p=tran+1;p<=end;p++)

    {

        i=p->i;

        j=p->j;

        r=p->a;

        x=p->b;

        t=p->c;

        b=t*(r*r+x*x);

        r/=b;

        x/=-b;

        b=t-1.0;

        ri=r*b;

        xi=x*b;

        rj=-ri/t;

        xj=-xi/t;

        vi=Um[i];

        vj=Um[j];

        vij=vi*vj;

        vi*=vi;

        vj*=vj;

        dij=Ua[i]-Ua[j];

        cd=vij*cos(dij);

        sd=vij*sin(dij);

        pij=vi*(ri+r)-r*cd-x*sd;

        pji=vj*(rj+r)-r*cd+x*sd;

        dpb=pij+pji;

        ph+=dpb;

        qij=-vi*(xi+x)+x*cd-r*sd;

        qji=-vj*(xj+x)+x*cd+r*sd;

        dqb=qij+qji; qh+=dqb;

        fprintf(fou," %3d  %3d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f\\n

            i,j,pij,qij,pji,qji,dpb,dqb);

    }

    fprintf(fou,"\\n\\n     系统总损耗:-   有功功率:%8.5f    无功功率:%8.5f

        ph,qh);

}

// 分配二维双精度型数组

double **newSpaceDouble2(int n1,int n2)

{

    int i;

    double **a=new double*[n1+1];

for(i=0;i<=n1;i++)

        a[i]=new double[n2+1];

    return a;

}

// 释放二维双精度型数组

void deleteSpaceDouble2(double **a,int n1)

{

for(int i=0;i<=n1;i++)

        delete[] a[i];

    delete[] a;

}    

三、输出潮流结果

程序运行后在“潮流输出.txt”文件中输出迭代过程的信息以及节点和支路潮流。

对于[例4-3]中的电力系统，输出结果如下：

      迭代次数:  1        最大功率误差:   0.500000

      迭代次数:  2        最大功率误差:   0.052346

      迭代次数:  3        最大功率误差:   0.001633

      迭代次数:  4        最大功率误差:   0.000002

                     *-*-*- 潮 流 计 算 结 果 *-*-*-

                        ——— 节 点 潮 流 ———

    no.i      Um         Ua         PG         QG         PL         QL

      1    0.96950   -3.87378    0.00000    0.00000    0.00000    0.00000

      2    1.03877   -9.23044    0.00000    0.00000    0.50000    0.30000

      3    1.05000   -1.84120    0.20000    0.19687    0.00000    0.00000

      4    1.05000    0.00000    0.47769    0.14430    0.15000    0.10000

                        ——— 支 路 潮 流 ———

   i    j       Pij        Qij        Pji        Qji        dP         dQ

   4    3    0.057   -0.05702   -0.05553    0.00179    0.00094   -0.05523

   1    4   -0.25735   -0.11014    0.27121    0.10132    0.01386   -0.00882

   2    2    0.00000   -0.05395    0.00000    0.00000    0.00000   -0.05395

   3    1    0.25553    0.19508   -0.24265   -0.19696    0.01288   -0.00187

   1    2    0.50000    0.30710   -0.50000   -0.24605    0.00000    0.06105

系统总损耗:-   有功功率: 0.02769    无功功率:-0.05883

另外，为了调试程序方便起见，在“中间结果.txt”文件中还输出了“导纳矩阵中的非零元素”，以及前两次迭代过程中的“雅可比矩阵中的非零元素”、“电压相位和有效值修正量”和“电压相位和有效值新值”等中间结果，可以用它们与[例4-3]中所给出的中间结果进行对比，以便核对或查找程序中存在的错误所在。

读者如有兴趣，还可以进一步作以下工作：

1.将输入数据中所有的节点编号改成节点名称(汉字和/或汉语拼音)，让程序自动进行节点编号；

2.加入节点编号优化；

3.对节点导纳矩阵和雅可比矩阵采用稀疏技巧存储和运算；

4.改成直角坐标形式的牛顿-拉夫逊法潮流计算程序。