潮流计算C语言编程-9节点 华北电力大学

辽宁工程技术大学电力系统分析课程设计设计题目9节点电力网络潮流计算指导教师院（系、部）专业班级学号姓名日期电气工程系课程设计标准评分模板电力系统分析课程设计任务书9节点系统单线图如下：基本数据如下：母线名基准电压区域号电压上限电压下限发电 1 16.5000 2 18.1500 14.8500 发电 2 18.000 1 19.800 16.2000 发电 3 13.8000 1 15.1800 12.4200 GEN1-230 230.000 2 0.0000 0.0000 GEN2-230 230.000 1 0.0000 0.0000 GEN3-230 230.000 1 0.0000 0.0000 STNA-230 230.000 2 0.0000 0.0000 STNB-230 230.000 2 0.0000 0.0000 STNC-230 230.000 1 0.0000 0.0000数据组I 侧母线J 侧母线编号所属区域单位正序电阻正序电抗正序充电电纳的1/2常规GEN1-230 STNA-230 1 I侧标么0.010000 0.085000 0.088000 常规STNA-230 GEN2-230 2 I侧标么0.032000 0.161000 0.153000 常规GEN2-230 STNC-230 3 I侧标么0.008500 0.072000 0.074500 常规STNC-230 GEN3-230 4 I侧标么0.011900 0.100800 0.104500 常规GEN3-230 STNB-230 5 I侧标么0.039000 0.170000 0.179000 常规STNB-230 GEN1-230 6 I侧标么0.017000 0.092000 0.079000表3 两绕组变压器数据负荷数据电网12-1班数据目录1 PSASP软件简介 (1)1.1 PSASP平台的主要功能和特点 (7)1.2 PSASP的平台组成 (8)2 牛顿拉夫逊潮流计算简介 (9)2.1 牛顿—拉夫逊法概要 (9)2.2 直角坐标下的牛顿—拉夫逊潮流计算 (11)2.3 牛顿—拉夫逊潮流计算的方法 (6)3 九节点系统单线图及元件数据 (8)3.1 九节点系统单线图 (8)3.2 系统各项元件的数据 (9)4 潮流计算的结果 (11)4.1 潮流计算后的单线图 (17)4.2 潮流计算结果输出表格 (18)5 结论 (22)6 参考文献 (17)1 PSASP软件简介“电力系统分析综合程序”（Power System Analysis Software Package,PSASP）是一套历史悠久、功能强大、使用方便的电力系统分析程序，是高度集成和开发具有我国自主知识产权的大型软件包。

(完整)C语言进行潮流计算编辑整理：尊敬的读者朋友们：这里是精品文档编辑中心，本文档内容是由我和我的同事精心编辑整理后发布的，发布之前我们对文中内容进行仔细校对，但是难免会有疏漏的地方，但是任然希望（(完整)C语言进行潮流计算）的内容能够给您的工作和学习带来便利。

同时也真诚的希望收到您的建议和反馈，这将是我们进步的源泉，前进的动力。

本文可编辑可修改，如果觉得对您有帮助请收藏以便随时查阅，最后祝您生活愉快业绩进步，以下为(完整)C语言进行潮流计算的全部内容。

电力系统课程设计C语言潮流计算学院：电气工程班级：电092班学号:0912002020学生姓名: 闵凯2013.3.7电力系统的潮流计算是对电力系统分析的最基本步骤也是最重要的步骤,是指在一定的系统结构和运行条件下，确定系统运行状态的计算,也即是对各母线（节点）电压,各元件（支路)传输电线或功率的计算。

通过计算出的节点电压和功率分布用以检查系统各元件是否过负荷，各点电压是否合理，以及功率损耗等。

即使对于一个简单的电力系统,潮流计算也不是一件简单就可以完成的事,其运算量很大,因此如果对于一个大的、复杂的电网来说的话，由于其节点多,分支杂，其计算量可想而知，人工对其计算也更是难上加难了。

特别是在现实生活中,遇到一个电力系统不会像我们期望的那样可以知道它的首端电压和首端功率或者是末端电压和末端功率，而是只知道它的首端电压和末端功率，更是使计算变的头疼万分。

为了使计算变的简单，我们就可以利用计算机，用C 语言编程来实现牛顿—拉夫逊（Newton —Raphson ）迭代法，最终实现对电力系统潮流的计算。

一．用牛顿—拉夫逊迭代法进行电力系统潮流计算的相关概念1.节点导纳矩阵如图所示的电力网络,将节点i 和j 的电压用•U i 和•.U j 表示,它们之间的支路导纳表示为y ij ,那么有基尔霍夫电流定律可知注入接点I 的电流•.I i （设流入节点的电流为正）等于离开节点I 的电流之和，因此有•i I)(.00••≠=≠=••-==∑∑j i nij ijnij iji U U II y （1-1）∴ •≠=≠=••∑∑-=U y y U I nij ij n ij ij i 00 （1-2)如令ii nij ij Y y =∑≠=0 ij ij Y y =-则可将（1-2）改写为:∑≠=••=nij ij ij i U Y I 1 I=1，2，…,n. （1-3)上式也可以写为: I =YU (1-4）其中Y 为节点导纳矩阵,也称为稀疏的对称矩阵,它是n ×n 阶方阵。

如图所示系统，试计算潮流分布，相关数据见《版潮流计算用户手册》P121。

#include<>#include <>float divRe(float b1,float b2,float b3,float b4){float a1r;a1r=(b1*b3+b2*b4)/(b3*b3+b4*b4);return(a1r);}float divIm(float b1,float b2,float b3,float b4){float a1i;a1i=(b2*b3-b1*b4)/(b3*b3+b4*b4);return(a1i);}float mulRe(float b1,float b2,float b3,float b4){float a2r;a2r=b1*b3-b2*b4;return(a2r);}float mulIm(float b1,float b2,float b3,float b4){float a2i;a2i=b2*b3+b1*b4;return(a2i);}float Max(float a[],int n){int i;float max;max= fabs(a[0]);for(i=1;i<n;i++){if(fabs(a[i])>max)max=fabs(a[i]);}return(max);}void main(){int i,j,k,h,km;int T=16;float eps,sumpi1,sumpi2,sumqi1,sumqi2,max,sumir,sumii,I1r,I1i,t,xx,xxx; float pi0[8],qi0[8],detpi[8],detqi[8],Iir0[8],Iii0[8],J0[16][16], detsi[16],detui[16], delta_p[9][9],delta_q[9][9],a[16][32],ni[16][16],H[8][8],N[8][8],J[8][8],L[8][8],ei1[9],fi1[9],sp[9][9],sq[9][9];static float ybr[9][9]={{,,0,0,0,,0,0,0},{,,,0,0,0,0,0,0},{0,,,,0,0,0,0,0},{0,0,,,,0,0,0,0},{0,0,0,,,,0,0,0},{,0,0,0,,,0,0,0},{0},{0},{0}};static float ybi[9][9]={{,,0,0,0,,0,0,},{,,,0,0,0,0,0,0},{0,,,,0,0,16,0,0},{0,0,,,,0,0,0,0},{0,0,0,,,,0,,0},{,0,0,0,,,0,0,0},{0,0,16,0,0,0,-16,0,0},{0,0,0,0,,0,0,,0},{,0,0,0,0,0,0,0,}};static float yd[9][9]={{0,,0,0,0,,0,0,0},{,0,,0,0,0,0,0,0},{0,,0,,0,0,0,0,0},{0,0,,0,,0,0,0,0},{0,0,0,,0,,0,0,0},{,0,0,0,,0,0,0,0},{0},{0},{0}} ;float ei0[9]={,,,,,,,,};float fi0[9]={};float pi[9]={0,,0,,0,,,,0};float qi[9]={0,,0,,0,,0,0,0};h=0;km=15;eps=;do{h+=1;printf("\nNow The %dth times\n",h);for(i=0;i<8;i++){printf("ei0[%d]=%f\t",i,ei0[i]);printf("fi0[%d]=%f\t",i,fi0[i]);}for(i=0;i<8;i++){printf("pi[%d]=%f\t",i,pi[i]);printf("qi[%d]=%f\t",i,qi[i]);}sumpi2=0;sumqi2=0;for(i=0;i<8;i++){for(j=0;j<9;j++){sumpi1=(ei0[i]*(ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j])+fi0[i]*(ybr[i][j]*fi0[j]+ybi [i][j]*ei0[j]));sumpi2+=sumpi1;}pi0[i]=sumpi2;printf("pi0[%d]=%f\t",i,pi0[i]);sumpi2=0;}for(i=0;i<8;i++){for(j=0;j<9;j++){sumqi1=(fi0[i]*(ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j])-ei0[i]*(ybr[i][j]*fi0[j]+ybi [i][j]*ei0[j]));sumqi2+=sumqi1;}qi0[i]=sumqi2;printf("qi0[%d]=%f\t",i,qi0[i]);sumqi2=0;}for(i=0;i<8;i++){detpi[i]=pi[i]-pi0[i];detqi[i]=qi[i]-qi0[i];if(i==6||i==7){qi0[i]=ei0[i]*ei0[i]+fi0[i]*fi0[i];detqi[i]=[i];}printf("detpi[%d]=%f\t",i,detpi[i]);printf("detqi[%d]=%f\t",i,detqi[i]);}//***************************************************************************** ********************//节点的注入电流表达式for(i=0;i<8;i++){ Iii0[i]=0;Iir0[i]=0;}for(i=0;i<8;i++){for(j=0;j<9;j++){Iir0[i]+=ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j];Iii0[i]+=ybr[i][j]*fi0[j]+ybi[i][j]*ei0[j];}}//***************************************************************************** ****************//求解NHJL矩阵for(i=0;i<8;i++){for(j=0;j<8;j++){if(i==j){if(i==6||i==7){H[i][j]=-ybi[i][j]*ei0[j]+ybr[i][j]*fi0[j]+Iii0[i];N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j]+Iir0[i];J[i][j]=2*fi0[i];L[i][j]=2*ei0[i];}else{H[i][j]=-ybi[i][j]*ei0[j]+ybr[i][j]*fi0[j]+Iii0[i];N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j]+Iir0[i];J[i][j]=-ybr[i][j]*ei0[j]-ybi[i][i]*fi0[j]+Iir0[i];L[i][j]=-ybi[i][j]*ei0[j]+ybr[i][j]*fi0[j]-Iii0[i];}}else{if(i==6||i==7){H[i][j]=ybr[i][j]*fi0[j]-ybi[i][j]*ei0[j];N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j];J[i][j]=0;L[i][j]=0;}else{ H[i][j]=ybr[i][j]*fi0[j]-ybi[i][j]*ei0[j];N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j];J[i][j]=-ybi[i][j]*fi0[j]-ybr[i][j]*ei0[j];L[i][j]=ybr[i][j]*fi0[j]-ybi[i][j]*ei0[j];}}}}//形成jacobian矩阵for(i=0;i<16;i++)for(j=0;j<16;j++){if(i%2==0&&j%2==0) J0[i][j]=H[i/2][j/2];else if(i%2==0&&j%2!=0) J0[i][j]=N[i/2][(j-1)/2];else if(i%2!=0&&j%2==0) J0[i][j]=J[(i-1)/2][j/2];else J0[i][j]=L[(i-1)/2][(j-1)/2];}//for(i=0;i<16;i++)//for(j=0;j<16;j++)// printf("J0[%d][%d]=%f\t",i,j,J0[i][j]);//***************************************************************************** *********//求detui//***************************************************************************** ********for(i=0;i<16;i++){if(i%2==0) detsi[i]=detpi[i/2];else detsi[i]=detqi[(i-1)/2];}//将detp和detq用一个数组表示for(i=0;i<T;i++)for(j=0;j<(2*T);j++){ if (j<16)a[i][j]=J0[i][j];else if (j==T+i)a[i][j]=;elsea[i][j]=;}for(i=0;i<T;i++){for(k=0;k<T;k++){if(k!=i){t=a[k][i]/a[i][i];for(j=0;j<(2*T);j++){xx=a[i][j]*t;a[k][j]=a[k][j]-xx;}}}}for(i=0;i<T;i++){t=a[i][i];for(j=0;j<(2*T);j++)a[i][j]=a[i][j]/t;}for(i=0;i<T;i++)for(j=0;j<T;j++)ni[i][j]=a[i][j+T];/*printf("逆矩阵为：\n");for (i=0;i<T;i++){for (j=0;j<T;j++)printf("%10.3f",ni[i][j]);printf("\n");} */xxx=;for(i=0;i<T;i++){xxx=;for(j=0;j<T;j++){xxx=xxx+ni[i][j]*detsi[j];}detui[i]=xxx;}//检测detui满足要求与否max=Max(detui,16);printf("max=%f\n",max);//****************************************************************** //算第n次迭代后的ufor(i=0;i<T;i++){if(i%2==0){ fi1[i/2]=fi0[i/2]+detui[i];}else{ei1[(i-1)/2]=ei0[(i-1)/2]+detui[i];}}//*********************************************for(i=0;i<8;i++)//下一次迭代赋初值{ei0[i]=ei1[i];fi0[i]=fi1[i];}for(i=0;i<8;i++){printf("ei0[%d]=%f\t",i,ei0[i]);printf("fi0[%d]=%f\t",i,fi0[i]);}for(i=0;i<8;i++){pi[i]=detpi[i]+pi0[i];qi[i]=detqi[i]+qi0[i];}}while(max>eps&&h<km);printf("All do %d times\n",h);sumir=0;sumii=0;//**********************************************for(i=0;i<9;i++)//平衡节点功率计算{I1r=mulRe(ybr[8][i],-ybi[8][i],ei0[i],-fi0[i]); I1i=mulIm(ybr[8][i],-ybi[8][i],ei0[i],-fi0[i]); sumir+=I1r;sumii+=I1i;}pi[8]=mulRe(ei0[8],fi0[8],sumir,sumii);qi[8]=mulIm(ei0[8],fi0[8],sumir,sumii);printf("S9=%f+j%f\n",pi[8],qi[8]);sumpi1=0;sumpi2=0;sumqi1=0;sumqi2=0;for(i=0;i<9;i++){for(j=0;j<9;j++)if(i!=j&&ybi[i][j]!=0){sumpi1=mulRe(ei0[i],fi0[i],,yd[i][j]);sumqi1=mulIm(ei0[i],fi0[i],,yd[i][j]);sumpi2=mulRe(ei0[i]-ei0[j],fi0[i]-fi0[j],-ybr[i][j],-ybi[i][j]);sumqi2=mulIm(ei0[i]-ei0[j],fi0[i]-fi0[j],-ybr[i][j],-ybi[i][j]);sumpi1+=sumpi2;sumqi1+=sumqi2;sp[i][j]=mulRe(ei0[i],fi0[i],sumpi1,-sumqi1);sq[i][j]=mulIm(ei0[i],fi0[i],sumpi1,-sumqi1);printf("S%d=%f+j%f\n",(i+1)*10+j+1,sp[i][j],sq[i][j]);}}for(i=0;i<9;i++){for(j=0;j<9;j++)if(j !=i&&ybi[i][j]!=0){delta_p[i][j]=sp[i][j]+sp[j][i];delta_q[i][j]=sq[i][j]+sq[j][i];if( (delta_p[i][j] !=0) || (delta_q[i][j] !=0)){printf("dS%d=%f+j%f\n",(i+1)*10+j+1,delta_p[i][j],delta_q[i][j]);}}}ei1[8]=ei0[8];fi1[8]=fi0[8];for(i=0;i<9;i++){printf("u%d=%f<%f\n",i+1,sqrt(ei1[i]*ei1[i]+fi1[i]*fi1[i]),atan(fi1[i]/ei1[i])*180/;} }。

//////////////////////////////////////////////////////////////////////// PQ分解法潮流////文件输入格式：节点总数n(包括联络节点),支路数zls ////节点数(发电机和负荷)nb,接地电抗数mdk,迭代精度eps // //考虑负荷静特性标志kk2(0考虑),平衡节点号，优化标志(0不优化) ////最大迭代次数it1,支路左右节点号izl[],jzl[],支路电阻zr[],电抗zx[] ////支路容纳zyk[],节点号nob[]及标志nobt[](0-PQ,-1-PV) ////发电机和负荷有功、无功pg[],qg[],pl[],ql[] ////电压v0[](pv节点输入实际值，PQ节点任输入一值) // //电抗节点号idk[],电抗值dkk[] ////////////////////////////////////////////////////////////////////////#include "math.h"#include "stdio.h"#define NS 2000 //最大节点数#define NS2 NS * 2#define NS4 1000 //NS4、NS必须大于2*zls。

#define ZS 3000 //最大支路数#define ZS2 ZS * 2#define DKS 200 //最大电抗器数#define N2 ZS * 4#define N3 ZS * 8 + NS * 4FILE *fp1, *fp2;char inname[12], outname[12];// fp1输入数据文件指针fp2输出文件指针// inname[]输入数据文件名outname[]输出数据文件名int n, zls, nb, mdk, mpj, bnsopton, it1, dsd, kk2, nzls;// 节点总数n(包括联络节点) 支路数(回路数)zls 节点数nb(发电机和负荷) // 接地电抗数mdk 精度eps 平衡节点号mpj// 节点优化(标志)bnsopton(=0节点不优化，!=0节点优化)// 最大迭代次数it1 最低电压或最大功率误差节点号dsd// 负荷静特性标志(=0考虑负荷静特性)// 支路数(双回线算一条支路)int izl[ZS], jzl[ZS], idk[DKS], yds[NS], ydz[NS], iy[ZS2];// izl[],jzl[],idk[]:分别存放左、右节点号和电抗器节点号。

C语言进行潮流计算C语言是一种通用的编程语言，被广泛用于开发各种应用程序。

虽然它的历史可以追溯到上个世纪70年代，但C语言在当今仍然非常流行。

它被广泛用于系统编程、嵌入式系统开发、游戏开发和科学计算等领域。

本文将重点介绍C语言的潮流计算。

潮流计算是一种用于分析电力系统稳态和动态行为的工具和技术。

它可以用于预测电力系统中潮流（电压和电流）的分布和变化。

潮流计算的目标是找到系统中各个节点的电压和相应的功率流向。

这对于电力系统的规划和运行至关重要。

C语言可以通过编写相应的程序来实现潮流计算。

在这个程序中，需要定义一组节点和支路，以及相应的导纳矩阵（节点之间的电导和电纳）。

导纳矩阵描述了电力系统中节点之间的电压和电流关系。

一旦定义了这些参数，就可以使用C语言编写算法来解决潮流计算问题。

潮流计算通常使用迭代的方法来求解。

在每一次迭代中，系统的节点电压和功率分布都会进行更新，直到达到稳定状态。

C语言提供了丰富的控制结构和数据类型，使得编写这样的迭代算法变得相对容易。

例如，可以使用for循环来实现迭代过程，同时使用浮点数类型来存储节点电压和功率分布。

潮流计算的核心是求解节点电压的方程组。

这个方程组通常是一个非线性方程组，需要使用牛顿-拉夫逊方法或高斯-赛德尔方法等迭代算法进行求解。

C语言提供了一些数值计算库，如GSL（GNU科学库），可以用于实现这些迭代算法。

这些库提供了各种数值计算函数，如求解非线性方程组、矩阵运算和插值等，可以大大简化潮流计算程序的编写。

除了求解节点电压方程组，潮流计算还需要考虑各种系统参数和条件。

例如，支路的参数（如电阻和电感）、发电机的容量和调节策略、负荷的需求和变化等。

C语言可以通过定义和存储这些参数，并在计算过程中引用它们来实现对这些条件的考虑。

此外，潮流计算还需要处理一些特殊情况，如无源系统和有环路系统。

无源系统是指没有外部电压源的电力系统，只由发电机和负荷组成。

有环路系统是指电力系统中存在环路，这会导致方程组无解。

《电力系统潮流上机》课程设计报告院系：电气与电子工程学院班级：电气1006学号学生姓名：指导教师：***设计周数：两周成绩：日期：2013 年1 月10 日一、 目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、 主要内容1． 手算.要求应用牛顿-拉夫逊法或P-Q 分解法手算求解，要求精度为0.001MW 。

节点1为平衡节点，电压︒∠=00.11U ，节点2为PQ 节点，负荷功率6.08.0~2j S +=，节点3是PV 节点，1.1,4.033==U P ，两条支路分别为04.001.013j Z +=，2.005.012j Z +=，对地支路33.030j y =。

2． 计算机计算(详见附录)3． 思考题3.1潮流计算的方法有哪些？各有何特点？答：潮流计算分为手算和机算两大类，其中机算又有高斯-赛德尔迭代法、牛顿-拉夫逊迭代法、P-Q 分解法等算法。

迭特点：手算求解潮流一般只用在简单的网络中，其计算量大，对于多节点的网络用手算一般难以解决问题，但通过手算可以加深物理概念的理解，还可以在运用计算机计算前以手算求取某些原始数据。

高斯-赛德尔迭代法：算法简单，对初值的要求不高，但需要迭代的次数多，收敛的速度慢，在早期的潮流计算程序中应用很多，之后逐渐被牛顿-拉夫逊迭代法所取代，但仍可作为计算程序前几次迭代的算法，以弥补后者对初值要求高的缺点。

牛顿-拉夫逊迭代法：是常用的解非线性方程组的方法，也是当前广泛采用的计算潮流的方法，其收敛速度快，几次迭代就可以得到最终的结果。

但其缺点是要求初值的选择得比较接近它们的精确值，否则迭代过程可能不收敛。

P-Q 分解法潮流计算：派生于以极坐标表示时的牛顿-拉夫逊法，其根据电力系统的特点，对后者的修正方程做了简化，P-Q 分解法的系数矩阵B ’和B ”代替了牛拉法中的雅可比矩阵J ，阶数降低，其中的元素在迭代过程中不发生变化，而且元素对称，这些都大大提高了运算速度，而且精确度几乎不受影响。

＃include <iostream〉#include 〈fstream>＃include<iomanip>#include〈math。

h〉using namespace std；//节点号类型负荷有功负荷无功母线数据（类型1＝PV节点，2＝PQ节点，3＝平衡节点）struct BUS｛int busno；int type；float Pd;float Qd;｝;//发电机数据节点号有功发电电压幅值struct Generator{int busno;float Pg;float Vg；｝;//支路信息节点I 节点J R X B/2 kstruct Line{int busi;int busj;float R;float X；float B;float k；｝；//deltaP deltaQ deltaV^2//void fun1（double YG[][50],double YB［][50]，double e［］，double f［］,int type[］，int N，double W［]，double P［],double Q［］，double V[］）｛double dP=0，dQ=0,dV=0;int i,j;for（i=0；i<N-1；i++)｛double A=0,B=0；for(j=0；j〈N；j++）{A+=YG［i]［j]＊e[j］-YB［i］[j］*f［j];B+=YG［i］［j］＊f[j]+YB［i][j]＊e[j］;}dV=V[i]＊V[i］-e［i]*e［i］—f［i］*f［i];dP=P[i]—e［i］＊A—f［i］＊B；W［2＊i]=dP;dQ=Q[i］-f［i]*A+e［i］＊B;if(type[i］==1)W[2＊i+1］=dQ;else W[2＊i+1]=dV；｝｝//Jacobi矩阵//void Jacobi（double YG［]［50］，double YB[］[50]，double e[50］,double f［50],int type［50]，int N ，double Ja［100］[101]）｛int i，j；for（i=0;i<N;i++)｛for(j=0；j〈N；j++){if(i！=j）{if（type［i]==1){Ja[2*i][2＊j]=—(YG[i]［j］＊e［i］+YB［i][j］＊f［i]）;Ja［2＊i]［2＊j+1]=YB［i]［j]＊e［i]-YG［i］［j］*f[i］;Ja［2＊i+1]［2*j]=Ja［2＊i］［2＊j+1］;Ja［2*i+1］[2*j+1]=—Ja［2*i］[2＊j]；}else {Ja［2*i］[2*j］=—YG［i］［j］*e[i］+YB[i][j］*f［i］;Ja［2＊i][2*j+1]=YB[i］[j］*e[i］-YG［i]［j]＊f［i]；Ja[2＊i+1]［2*j+1]=Ja[2＊i+1］［2＊j]=0；}｝else {double a[50］={0｝，b[50]=｛0};for（int k=0；k<N；k++){a[i]+=(YG[i］［k]*e[k]—YB［i]［k］*f[k］）;b[i]+=（YG[i][k]*f[k]+YB[i］[k]*e［k］）；Ja［2＊i]［2*j］=—a［i]-YG［i]［i］*e［i］-YB[i]［i]＊f［i］；Ja[2*i]［2*j+1］=-b[i］+YB[i]［i］*e[i]-YG［i］［i]＊f[i]；if(type［i]==1）{Ja[2*i+1][2＊j］=b［i]+YB[i］［i]＊e［i］—YG［i］[i]＊f[i]；Ja[2＊i+1][2*j+1]=-a［i]+YG［i]［i]＊e［i]+YB［i］[i］＊f[i];}else ｛Ja[2*i+1］［2*j］=—2＊e[i]；Ja［2*i+1]［2*j+1］=-2＊f[i];}}｝｝}}//高斯消元法解方程组函数//void gauss（double a[］［101］，int n）｛int i，j,k；double c;for（k=0；k〈n-1；k++) ｛c=a[k］[k］；for（j=k；j〈=n；j++) a［k][j］/=c；for（i=k+1；i〈n;i++）｛c=a［i][k］；for(j=k；j<=n；j++）a[i][j］—=c*a［k］[j];｝}a［n-1]［n]/=a[n-1][n-1];for（k=n-2;k>=0;k——）for(j=k+1；j<n；j++) a［k]［n]-=a[k］［j]＊a［j]［n］;}void main（){ifstream fin；int N=0，GS=0，LD=0,ZLs=0; //节点数发电机数负荷数支路数// BUS *B；Generator *G；Line ＊L;//从文本中读入原始数据到数组中//fin。

C 语言程序设计潮流计算学院自动化学院专业班级学号姓名联系方式本程序潮流计算部分采用牛顿拉夫逊极坐标法进行计算，求解一次多元方程采用高斯列主元分解法进行求解。

根据工程实际，在存储文件实时记录产生文件时间。

此外本程序特意增加了文件查看功能，方便文件的查看。

程序代码：#include<stdio.h>#include<malloc.h>#include<stdlib.h>#include<math.h>#include<string.h>#include<time.h>#include<conio.h>void dataprepare(void);void initial(void);void Yc(void);void showy(void);double detpqc(void);void showdetav(void);void Jrc(void);void showdetpq(void);void Gauss(void);void showj(void);void showsolution(void);char save2file(void);char option(void);int chose(void);int list(char filename[],int i,char dstn[]);struct PQV{char name[5];double vb;double p;double q;double v;double ag;}*pqv;struct Y{double G;double B;}**y;struct LZ{char name1[5];double vb1;char name2[5];double vb2;double r;double x;double b;}*lz;struct TZ{char name1[5];double vb1;char name2[5];double vb2;double x;double k1;double k2;}*tz;char fsource[20];int pqnum=0,pvnum=0,lznum=0,tznum=0;int temp1=0,temp2=0,temp3=0,temp4=0;int diedai=0;int numw=1,numofw=0;double *detpq,*detav,*Pi,*Qi,error,wucha,**J; /* detpq节点功率的误差量，detav修正量，Pi 节点的有功功率*/char c,ch,dstname[20],filename[20];int main(){opt:c=option();if(c=='1')goto opt1;else if(c=='2'){list("1.bin",1,filename);system("mode con:cols=100 lines=30 & color 07");while(c!='0'){system("cls");printf(" ################################计算结果文件查看################################\n");printf(" [8]上一个[2]下一个[5]查看[-]删除文件【0】退出\n\n");list("1.bin",0,filename);chose();goto opt;}else if(c=='0')exit(0);else goto opt;opt1:dataprepare();init:initial();Yc();printf("输入允许误差量：");scanf("%lf",&wucha);printf("选择模式：a.详情；b.简单：");redo:fflush(stdin);c=getchar();if(c!='a'&&c!='b'&&c!='A'&&c!='B'){printf("输入错误！重新输入：");goto redo;}if(c!='B'&&c!='b')showy();error=detpqc();showdetpq();while(error>wucha&&error<100){if(c!='B'&&c!='b'){printf("\t\t是否继续迭代<空格：是R/r:重赋初值E/e:退出潮流计算>？\n");ifdie: fflush(stdin);ch=getch();switch(ch){case ' ':break;case 'R':;case 'r':goto init;break;case 'E':;case 'e':goto opt;break;default: goto ifdie;}}diedai++;Jrc();if(c!='B'&&c!='b'){printf("\n>>>>>>>>>>>>>>>>>>>第%d 次迭代<<<<<<<<<<<<<<<<<<<<<\n",diedai);showj();Gauss();if(c!='B'&&c!='b')showdetav();for(temp1=0;temp1<pqnum+pvnum;temp1++) /* 修正电压*/pqv[temp1].ag=pqv[temp1].ag-detav[temp1];for(temp1=0;temp1<pqnum;temp1++)pqv[temp1].v=pqv[temp1].v-detav[temp1+pqnum+pvnum];error=detpqc();if(c!='B'&&c!='b')showdetpq();}if(error>=300||error<0){printf("\t\t迭代不收敛<R/r:重赋初值E/e:退出潮流计算>？\n");what: fflush(stdin);ch=getch();switch(ch){case 'R':;case 'r':goto init;break;case 'E':;case 'e':goto opt;break;default: goto what;}}showsolution();ch=save2file();if(ch=='1')goto opt;elsegoto init;}void dataprepare(void){FILE *fp;char ch,type[4];pqnum=0,pvnum=0,lznum=0,tznum=0,temp1=0,temp2=0,temp3=0,temp4=0;printf("\n 请输入潮流计算的数据文件\n --");getf:fflush(stdin);gets(fsource);while((fp=fopen(fsource,"r"))==NULL) /* 打开文件*/ {printf("\t\t不能打开文件！！\n\t\t重新输入请键入“Y”，其他任意键退出:");fflush(stdin);ch=getchar();if(ch=='Y'||ch=='y'){printf("\n\t\t潮流计算的数据文件：");goto getf;}elseexit(0);}system("mode con:cols=130 lines=50 & color 06");while(!feof(fp)) /* 计算各种节点数目*/{ch=fgetc(fp);if(ch==10){if(fgets(type,3,fp)){if(strcmp(type,"pq")==0)pqnum++;if(strcmp(type,"pv")==0)pvnum++;if(strcmp(type,"lz")==0)lznum++;if(strcmp(type,"tz")==0)tznum++;}}}pqv=(struct PQV*)malloc((pqnum+pvnum+1)*sizeof(struct PQV)); /* 根据节点数目开辟各节点储存空间*/lz=(struct LZ*)malloc((lznum)*sizeof(struct LZ));tz=(struct TZ*)malloc((tznum)*sizeof(struct TZ));rewind(fp);while(!feof(fp)) /* 读取和输入各节点数据*/{ch=fgetc(fp);if(ch==10){if(!fgets(type,3,fp))break;if(strcmp(type,"pq")==0){fscanf(fp," %s %lf %lf %lf",pqv[temp1].name,&pqv[temp1].vb,&pqv[temp1].p,&pqv[temp1]. q);temp1++;}if(strcmp(type,"pv")==0){fscanf(fp," %s %lf %lf\t%lf",pqv[temp2+pqnum].name,&pqv[temp2+pqnum].vb,&pqv[temp2 +pqnum].p,&pqv[temp2+pqnum].v);temp2++;}if(strcmp(type,"ph")==0){fscanf(fp," %s %lf\t\t%lf %lf",pqv[pqnum+pvnum].name,&pqv[pqnum+pvnum].vb,&pqv[pqn um+pvnum].v,&(pqv[pqnum+pvnum].ag));pqv[pqnum+pvnum].ag=pqv[pqnum+pvnum].ag/180*3.1415926;}if(strcmp(type,"lz")==0){fscanf(fp," %s %lf %s %lf %lf %lf %lf",lz[temp3].name1,&lz[temp3].vb1,lz[temp3].name2,&lz[ temp3].vb2,&lz[temp3].r,&lz[temp3].x,&lz[temp3].b);temp3++;}if(strcmp(type,"tz")==0){fscanf(fp," %s %lf %s %lf %lf %lf %lf",tz[temp4].name1,&tz[temp4].vb1,tz[temp4].name2,&tz [temp4].vb2,&tz[temp4].x,&tz[temp4].k1,&tz[temp4].k2);temp4++;}}}fclose(fp);puts("******************************************************************************* ********************************************");puts(" 电力系统潮流计算");puts("******************************************************************************* ********************************************");printf("从文件中获得的数据如下：\n");printf("节点类型名称额定电压有功标幺无功标幺电压标幺相角弧度\n");for(temp1=0;temp1<=pqnum+pvnum;temp1++){if(temp1<pqnum)printf("PQ %s %9.3lf %9.3lf %9.3lf\n",pqv[temp1].name,pqv[temp1].vb,pqv[temp1].p,pqv[temp1].q);else if(temp1<pqnum+pvnum)printf("PV %s %9.3lf %9.3lf %9.3lf\n",pqv[temp1].name,pqv[temp1].vb,pqv[te mp1].p,pqv[temp1].v);elseprintf(" 平衡%s %9.3lf %9.3lf %9.3lf\n",pqv[temp1].name,pqv[temp 1].vb,pqv[temp1].v,pqv[temp1].ag);}printf("\n");printf(" 线路节点1 额定电压节点2 额定电压电阻电抗对地导纳\n");for(temp1=0;temp1<lznum;temp1++)printf("输电线路%s %9.3lf %s %9.3lf %9.3lf %9.3lf %9.5lf\n",lz[temp1].name1,lz[temp1].vb1,lz[te mp1].name2,lz[temp1].vb2,lz[temp1].r,lz[temp1].x,lz[temp1].b);printf("\n");printf(" 线路节点1 额定电压节点2 额定电压电抗值一次比二次\n");for(temp1=0;temp1<tznum;temp1++)printf(" 变压器%s %9.3lf %s %9.3lf %9.3lf %9.3lf ：%9.3lf\n",tz[temp1].name1,tz[temp1].vb1, tz[temp1].name2,tz[temp1].vb2,tz[temp1].x,tz[temp1].k1,tz[temp1].k2);}void initial(void){char c;re_insert:for(temp1=0;temp1<pqnum+pvnum;temp1++){if(temp1<pqnum){printf("请输入第%d个节点%s(pq)的电压幅值：",temp1+1,pqv[temp1].name);fflush(stdin);scanf("%lf",&pqv[temp1].v);printf("请输入第%d个节点%s(pq)的电压相角：",temp1+1,pqv[temp1].name);fflush(stdin);scanf("%lf",&pqv[temp1].ag);pqv[temp1].ag=pqv[temp1].ag/180*3.1415926;}if(temp1>=pqnum){printf("请输入第%d个节点%s(pv)的电压相角：",temp1+1,pqv[temp1].name);fflush(stdin);scanf("%lf",&pqv[temp1].ag);}}printf("\t\t\t按ESC重新赋值，其他任意键继续！\n");c=getch();if(c==27)goto re_insert;}void Yc(void){int i, j,temp,flag=0;y=(struct Y **)malloc((pqnum+pvnum+1)*sizeof(struct Y *));for(temp1=0;temp1<pqnum+pvnum+1;temp1++)y[temp1]=(struct Y *)malloc((pqnum+pvnum+1)*sizeof(struct Y));for(i=0;i<=pqnum+pvnum;i++){for(j=0;j<=(pqnum+pvnum);j++){if(i!=j){flag=0;for(temp=0;temp<lznum;temp++) /* 当两节点间为线路时*/if((strcmp(pqv[i].name,lz[temp].name1)==0&&strcmp(pqv[j].name,lz[temp].name2)==0)||(s trcmp(pqv[i].name,lz[temp].name2)==0&&strcmp(pqv[j].name,lz[temp].name1)==0)){ flag++;y[i][j].G=-lz[temp].r/(pow(lz[temp].r,2)+pow(lz[temp].x,2));y[i][j].B=lz[temp].x/(pow(lz[temp].r,2)+pow(lz[temp].x,2));}for(temp=0;temp<tznum;temp++) /* 当两节点间为变压器*/if((strcmp(pqv[i].name,tz[temp].name1)==0&&strcmp(pqv[j].name,tz[temp].name2)==0)||( strcmp(pqv[i].name,tz[temp].name2)==0&&strcmp(pqv[j].name,tz[temp].name1)==0)){ flag++;y[i][j].G=0;y[i][j].B=1/tz[temp].x*tz[temp].k1/tz[temp].k2;}if(flag==0){y[i][j].G=0;y[i][j].B=0;}if(flag>1)printf("\n两节点间出现多条支路，本程序暂不支持!!\n");}}}for(i=0;i<=pqnum+pvnum;i++){y[i][i].G=0,y[i][i].B=0;for(temp=0;temp<lznum;temp++) /* 当与另一节点间为线路时*/if((strcmp(pqv[i].name,lz[temp].name1)==0)||(strcmp(pqv[i].name,lz[temp].name2)==0)) {y[i][i].G+=lz[temp].r/(pow(lz[temp].r,2)+pow(lz[temp].x,2));y[i][i].B+=-lz[temp].x/(pow(lz[temp].r,2)+pow(lz[temp].x,2))+lz[temp].b/2;}for(temp=0;temp<tznum;temp++) /* 当与另一节点间为变压器*/{if(strcmp(pqv[i].name,tz[temp].name1)==0){y[i][i].B+=-(1/tz[temp].x*tz[temp].k1/tz[temp].k2+1/tz[temp].x*(tz[temp].k2-tz[temp].k1)/tz[ temp].k2);}if(strcmp(pqv[i].name,tz[temp].name2)==0){y[i][i].B+=-(1/tz[temp].x*tz[temp].k1/tz[temp].k2+1/tz[temp].x*(tz[temp].k1-tz[temp].k2)*tz[ temp].k1/pow(tz[temp].k2,2));}}}}void showy(void){printf("\n");printf("导纳矩阵Y:\n");for(temp1=0;temp1<=pqnum+pvnum;temp1++){for(temp2=0;temp2<=pqnum+pvnum;temp2++)printf("%9.6lf+j%9.6lf ",y[temp1][temp2].G,y[temp1][temp2].B);printf("\n");}printf("\n");}double detpqc(void){double max=0;detpq=(double *)malloc((2*pqnum+pvnum)*sizeof(double)); /* 为detpq误差量*/ detav=(double *)malloc((2*pqnum+pvnum)*sizeof(double)); /* 修正量*/Pi=(double *)malloc((pqnum+pvnum)*sizeof(double)); /* */Qi=(double *)malloc((pqnum+pvnum)*sizeof(double));for(temp1=0;temp1<2*(pqnum+pvnum);temp1++){if(temp1<pqnum+pvnum) /* 求P（pq,pv)误差量pqnum+pvnum个*/{Pi[temp1]=0;for(temp2=0;temp2<=(pqnum+pvnum);temp2++)Pi[temp1]=Pi[temp1]+pqv[temp2].v*(y[temp1][temp2].G*cos(pqv[temp1].ag-pqv[temp2].ag )+y[temp1][temp2].B*sin(pqv[temp1].ag-pqv[temp2].ag));Pi[temp1]=pqv[temp1].v*Pi[temp1];detpq[temp1]=pqv[temp1].p-Pi[temp1];if(fabs(detpq[temp1])>max)max=fabs(detpq[temp1]);}if(temp1>=pqnum+pvnum) /* 求Q（pq）误差量pqnum个*/{Qi[temp1-pqnum-pvnum]=0;for(temp2=0;temp2<=(pqnum+pvnum);temp2++)Qi[temp1-pqnum-pvnum]=Qi[temp1-pqnum-pvnum]+pqv[temp2].v*(y[temp1-pqnum-pvnu m][temp2].G*sin(pqv[temp1-pqnum-pvnum].ag-pqv[temp2].ag)-y[temp1-pqnum-pvnum][temp2 ].B*cos(pqv[temp1-pqnum-pvnum].ag-pqv[temp2].ag));Qi[temp1-pqnum-pvnum]=pqv[temp1-(pqnum+pvnum)].v*Qi[temp1-pqnum-pvnum];detpq[temp1]=pqv[temp1-(pqnum+pvnum)].q-Qi[temp1-pqnum-pvnum];if(fabs(detpq[temp1])>max&&temp1!=2*pqnum+pvnum)max=fabs(detpq[temp1]);}}return max;}void showdetav(void){printf("解得:\n");for(temp1=0;temp1<2*pqnum+pvnum;temp1++)if(temp1<pqnum)printf("Δδ%d：%lf ",temp1,detav[temp1]);elseprintf("ΔV %d: %lf ",temp1-pqnum,detav[temp1]);printf("\n");}void Jrc(void){int i,j;double agij;J=(double **)malloc((2*pqnum+pvnum)*sizeof(double *)); /* 为雅可比矩阵分配内存空间*/for(temp1=0;temp1<2*pqnum+pvnum;temp1++)J[temp1]=(double *)malloc((2*pqnum+pvnum)*sizeof(double));for(i=0;i<pqnum+pvnum;i++){for(j=0;j<pqnum+pvnum;j++) /* Hij */{if(i!=j){agij=pqv[i].ag-pqv[j].ag;J[i][j]=-pqv[i].v*pqv[j].v*(y[i][j].G*sin(agij)-y[i][j].B*cos(agij));}elseJ[i][j]=pqv[i].v*pqv[i].v*y[i][i].B+Qi[i];}for(j=0;j<pqnum;j++) /* Nij */{if(i!=j){agij=pqv[i].ag-pqv[j].ag;J[i][j+pqnum+pvnum]=-pqv[i].v*(y[i][j].G*cos(agij)+y[i][j].B*sin(agij));}elseJ[i][j+pqnum+pvnum]=-pqv[i].v*y[i][i].G-Pi[i];}}for(i=0;i<pqnum;i++){for(j=0;j<pqnum+pvnum;j++) /* Kij */{if(i!=j){agij=pqv[i].ag-pqv[j].ag;J[i+pqnum+pvnum][j]=pqv[i].v*pqv[j].v*(y[i][j].G*cos(agij)+y[i][j].B*sin(agij));}elseJ[i+pqnum+pvnum][j]=pqv[i].v*pqv[i].v*y[i][i].G-Pi[i];}for(j=0;j<pqnum;j++) /* Lij */{if(i!=j){agij=pqv[i].ag-pqv[j].ag;J[i+pqnum+pvnum][j+pqnum+pvnum]=-pqv[i].v*(y[i][j].G*sin(agij)-y[i][j].B*cos(agij));}elseJ[i+pqnum+pvnum][j+pqnum+pvnum]=pqv[i].v*y[i][i].B-Qi[i];}}}void showj(void){printf("\n");printf("雅可比矩阵J:\n");for(temp1=0;temp1<2*pqnum+pvnum;temp1++){for(temp2=0;temp2<2*pqnum+pvnum;temp2++)printf("%9.6lf ",J[temp1][temp2]);printf("\n");}printf("\n");}void Gauss(void){int xiao,big,i,j,n=2*pqnum+pvnum;double biggest,sta,E;double conj=0,conpq=0;for(xiao=0;xiao<2*pqnum+pvnum-1;xiao++){big=xiao;biggest=fabs(J[xiao][xiao]);for(i=xiao+1;i<2*pvnum+pqnum;i++) /* 找出列主元*/ if(fabs(J[i][xiao])>biggest){big=i;biggest=fabs(J[i][xiao]);}for(i=0;i<2*pqnum+pvnum;i++) /* 交换两行*/{conj=J[xiao][i];J[xiao][i]=J[big][i];J[big][i]=conj;}conpq=detpq[xiao];detpq[xiao]=detpq[big];detpq[big]=conpq;for(i=xiao+1;i<n;i++) /* 消元*/{sta=J[i][xiao];for(j=0;j<n;j++)J[i][j]=J[i][j]-sta*J[xiao][j]/J[xiao][xiao];detpq[i]=detpq[i]-sta*detpq[xiao]/J[xiao][xiao];}}detav[n-1]=detpq[n-1]/J[n-1][n-1]; /* 回代*/for(i=n-2;i>=0;i--){ E=0;for(j=i+1;j<2*pqnum+pvnum;j++)E=E+J[i][j]*detav[j];detav[i]=(detpq[i]-E)/J[i][i];}}void showdetpq(void){printf("\n");printf("功率不平衡量:\n");for(temp1=0;temp1<2*pqnum+pvnum;temp1++){if(temp1<pqnum+pvnum)printf("ΔP%d: %lf ",temp1+1,detpq[temp1]);elsepr intf("ΔQ%d: %lf ",temp1-pqnum,detpq[temp1]);}printf("\n最大功率不平衡量为（绝对值）：%lf\n",error);printf("\n");}void showsolution(void){printf("\n>>>>>>>>>>>>>>>>>>> 迭代结束<<<<<<<<<<<<<<<<<<<<<\n");printf("\t\t总共进行%d次迭代:\n",diedai);for(temp1=0;temp1<pqnum+pvnum;temp1++)printf("节点%d<%s>电压：幅值%lf 相角%lf\n",temp1+1,pqv[temp1].name,pqv[temp1].v,pqv[temp1].ag);for(temp1=0;temp1<(pqnum+pvnum+1);temp1++){Pi[temp1]=0;for(temp2=0;temp2<=(pqnum+pvnum);temp2++){Pi[temp1]=Pi[temp1]+pqv[temp2].v*(y[temp1][temp2].G*cos(pqv[temp1].ag-pqv[temp2].ag )+y[temp1][temp2].B*sin(pqv[temp1].ag-pqv[temp2].ag));Qi[temp1]=Qi[temp1]+pqv[temp2].v*(y[temp1][temp2].G*sin(pqv[temp1].ag-pqv[temp2].a g)-y[temp1][temp2].B*cos(pqv[temp1].ag-pqv[temp2].ag));}Pi[temp1]=Pi[temp1]*pqv[temp1].v;Qi[temp1]=Qi[temp1]*pqv[temp1].v;printf("节点%d<%s>功率: %8.6lf+j%8.6lf\n",temp1+1,pqv[temp1].name,Pi[temp1],Qi[temp1]);}printf("\n");}char save2file(void){FILE *fp;int i;time_t t;char filename[20],c;printf("\t\t是否保存文件<Y/y>：");fflush(stdin);c=getchar();if(c!='Y'&&c!='y')goto next;printf("输入要储存数据的文件名(默认后缀.data)：");fflush(stdin);scanf("%s",filename);strcat(filename,".data");if((fp=fopen(filename,"w"))==NULL){printf("无法执行存储操作！！");exit(0);}time(&t);fprintf(fp,"\t\t\t潮流计算结果\n保存时间：%s\n对应数据文件：%s\n",ctime(&t),fsource);fprintf(fp,"<标幺值>\n");for(i=0;i<pqnum;i++)fprintf(fp," PQ节点%s，电压幅值：%10.6lf 相角:%10.6lf度有功功率：%10.6lf 无功功率：%10.6lf\n",pqv[i].name,pqv[i].v,pqv[i].ag*180/3.14159265337,Pi[i],Qi[i]);for(i=pqnum;i<pqnum+pvnum;i++)fprintf(fp," PV节点%s，电压幅值：%10.6lf 相角:%10.6lf度有功功率：%10.6lf 无功功率：%10.6lf\n",pqv[i].name,pqv[i].v,pqv[i].ag*180/3.14159265337,Pi[i],Qi[i]);fprintf(fp,"平衡节点%s，电压幅值：%10.6lf 相角:%10.6lf度有功功率：%10.6lf 无功功率：%10.6lf\n",pqv[i].name,pqv[i].v,pqv[i].ag*180/3.14159265337,Pi[i],Qi[i]);fclose(fp);printf("\t\t\t文件保存成功\n");next:puts("\t1.返回主界面 2.重新计算此文件数据:");nex:fflush(stdin);c=getch();if(c!='1'&&c!='2')goto nex;return c;}char option(void){char c;system("mode con:cols=75 lines=25 &color 16");cls:puts(" ***************************************************************");puts(" 电力系统潮流计算");puts("***************************************************************\n\n");printf(" 选择您的操作：\n\n\t\t1.潮流计算 2.查看数据文件0.退出程序\n\t\t\t\t");fflush(stdin);c=getch();return c;}int list(char filename[],int i,char dstn[]){FILE *fp;char ch,c;int shu=0,j=0,k=0;char name[50];system("dir /x *.data>1.bin 2>nul ");if((fp=fopen("1.bin","r"))==NULL){printf("读取列表有误！");system("pause");return 0;}while((ch=fgetc(fp))!=EOF){if(ch=='2'&&c==10){shu++;if(shu!=i){fseek(fp,47L,1);fgets(name,50L,fp);}if(shu==i){fseek(fp,35L,1);fgets(dstn,13,fp);if(dstn[0]==' '){fseek(fp,1L,1);fgets(dstn,13,fp);dstn[strlen(dstn)-1]=0;}return 0;}if(i==0){printf("\t\t");if(numw==shu)printf("\b\b->");printf(" %3d . %s",shu,name);numofw=shu;}ch=10;}c=ch;}fclose(fp);return 0;}int chose(void){char del[30]={"del "};FILE *fp;fflush(stdin);c=getch();if(c=='2'&&numw<numofw)numw++;if(c=='8'&&numw>1)numw--;if(c=='-'){list("1.bin",numw,filename);strcat(del,filename);strcat(del," 2>null");system(del);}if(c=='5') /* 5 is play key */{list("1.bin",numw,filename);if((fp=fopen(filename,"r"))==NULL){printf("读取列表有误！");system("pause");return 0;}while((ch=fgetc(fp))!=EOF)putchar(ch);system("pause & del /f /Q 2>nul");fclose(fp);}return 0;}运行效果：启动主界面：1.键入1，开始潮流计算，输入潮流计算数据文件：程序将自动读取数据，并进入如图界面，列出数据：此时需为潮流计算设定初值，如下图：若需重新给定初值则键入按ESC键，否则按任意键继续：按任意键给定允许误差，并选择显示模式：a.详情可列出计算过程的各种数据，如导纳矩阵，雅可比矩阵，误差量等。

0 概述电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各线的电压、各元件中流过的功率、系统的功率损耗等等。

在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性、可靠性和经济性。

此外，在进行电力系统静态及暂态稳定计算时，要利用潮流计算的结果作为其计算的基础；一些故障分析以及优化计算也需要有相应的潮流计算作配合；潮流计算往往成为上述计算程序的一个重要组成部分。

本设计首先简单介绍了电力系统分析综合程序——PSASP ,然后对牛顿—拉夫逊潮流计算法的原理进行了简单的讲解，再以一个9节点电力系统为例，介绍了PSASP在潮流计算中的应用，最后给出了利用PSASP计算潮流的结果。

1 PSASP—电力系统分析综合程序简介电力系统分析综合程序（Power System Analysis Software Package，PSASP）是中国电力科学研究院电力系统技术分公司（原电网数字仿真技术研究所）开发的一套用于进行电力系统分析计算的软件包，其主要包括如下模块：PSASP 图模一体化平台PSASP 潮流计算模块（LF）PSASP 暂态稳定计算模块（ST ）PSASP 短路计算模块（SC ）PSASP 最优潮流和无功优化计算模块（OPF ）PSASP 静态安全分析模块（SA ）PSASP 网损分析模块（NL）PSASP 静态和动态等值计算模块（EQ）PSASP 用户自定义模型和程序接口模块（UD/UPI）PSASP 直接法稳定计算模块（DST）PSASP 小干扰稳定分析模块（SST ）PSASP 电压稳定分析模块（VST ）PSASP 继电保护整定计算模块（RPS）PSASP 线性/非线性参数优化模块（LPO/NPO）PSASP 谐波分析模块（HMA）PSASP 分布式离线计算平台PSASP 电网风险评估系统PSASP 暂态稳定极限自动求解程序PSASP 负荷电流防冰融冰辅助决策系统PSASP 功能强大、使用方便、高度集成并开放，是具有我国自主知识产权的大型软件包。

// flow.cpp : Defines the entry point for the console application.//#include "stdafx.h"#include <iostream.h>#include <fstream.h>#include "NEquation.h"#include "math.h"#include "config.h"void test(){NEquation ob1;ob1.SetSize(2);ob1.Data(0,0)=1;ob1.Data(0,1)=2;ob1.Data(1,0)=2;ob1.Data(1,1)=1;ob1.V alue(0)=5;ob1.V alue(1)=4;ob1.Run();printf("x1=%f\n",ob1.V alue(0));printf("x2=%f\n",ob1.V alue(1));}void GetData()//Read the data{FILE *fp;int i;if((fp=fopen("F:\\1061180523\\flow\\data\\data.txt","r"))==NULL){printf("Can not open the file named 'data.txt' \n");return;}for(i=0;i<=Bus_Num-1;i++){fscanf(fp,"%d,%f,%f,%f,%f,%f,%f,%d",&gBus[i].No,&gBus[i].V oltage,&gBus[i].Phase, &gBus[i].GenP,&gBus[i].GenQ,&gBus[i].LoadP,&gBus[i].LoadQ,&gBus[i].Type);}for(i=0;i<=Line_Num-1;i++){fscanf(fp,"%d,%d,%d,%f,%f,%f,%f",&gLine[i].No,&gLine[i].No_I,&gLine[i].No_J, &gLine[i].R,&gLine[i].X,&gLine[i].B,&gLine[i].k);}fclose(fp);}void GetYMatrix(){int i,j;FILE *fp;// calculate the Y matrixint bus1,bus2;float r,x,d,g,b;for(i=0;i<=Bus_Num-1;i++) //导纳矩阵赋初值为零{for(j=0;j<=Bus_Num-1;j++){gY_G[i][j]=0;gY_B[i][j]=0;}}for(i=0; i<=Line_Num-1; i++){bus1=gLine[i].No_I-1; //线路一侧连接的节点号减1bus2=gLine[i].No_J-1; //线路另一侧连接的节点号减1r=gLine[i].R;x=gLine[i].X;d=r*r+x*x;g=r/d;b=-x/d;if(gLine[i].k==0) //没有变压器的情况{gY_G[bus1][bus1]=gY_G[bus1][bus1]+g;//gY_G[bus1][bus1]=ggY_G[bus2][bus2]=gY_G[bus2][bus2]+g;//gY_G[bus2][bus2]=ggY_G[bus1][bus2]=gY_G[bus1][bus2]-g;//gY_G[bus1][bus2]=-g 矩阵是对称的gY_G[bus2][bus1]=gY_G[bus2][bus1]-g;//gY_G[bus2][bus1]=-ggY_B[bus1][bus1]=gY_B[bus1][bus1]+b+gLine[i].B;//gY_B[bus1][bus1]=b+gLine[i].BgY_B[bus2][bus2]=gY_B[bus2][bus2]+b+gLine[i].B;//gY_B[bus2][bus2]=b+gLine[i].BgY_B[bus1][bus2]=gY_B[bus1][bus2]-b;//gY_B[bus1][bus2]=-bgY_B[bus2][bus1]=gY_B[bus2][bus1]-b;//gY_B[bus2][bus1]=-b}else //有变压器的情况{gY_G[bus1][bus1]=gY_G[bus1][bus1]+g/gLine[i].k+g*(1-gLine[i].k)/gLine[i].k/gLine[i].k;gY_G[bus2][bus2]=gY_G[bus2][bus2]+g/gLine[i].k+g*(gLine[i].k-1)/gLine[i].k;gY_G[bus1][bus2]=gY_G[bus1][bus2]-g/gLine[i].k;gY_G[bus2][bus1]=gY_G[bus2][bus1]-g/gLine[i].k;gY_B[bus1][bus1]=gY_B[bus1][bus1]+b/gLine[i].k+b*(1-gLine[i].k)/gLine[i].k/gLine[i].k;gY_B[bus2][bus2]=gY_B[bus2][bus2]+b/gLine[i].k+b*(gLine[i].k-1)/gLine[i].k;gY_B[bus1][bus2]=gY_B[bus1][bus2]-b/gLine[i].k;gY_B[bus2][bus1]=gY_B[bus2][bus1]-b/gLine[i].k;}}// output the Y matrixif((fp=fopen("F:\\1061180523\\flow\\data\\ymatrix.txt","w"))==NULL){printf("Can not open the file named 'ymatrix.txt' \n");return ;}fprintf(fp,"---Y Matrix---\n");for(i=0;i<=Bus_Num-1;i++){for(j=0;j<=Bus_Num-1;j++){fprintf(fp,"Y(%d,%d)=(%10.5f,%10.5f)\n",i+1,j+1,gY_G[i][j],gY_B[i][j]);}}fclose(fp);}void SetInitial(){int i;for(i=0;i<=Bus_Num-1;i++){if(gBus[i].Type==3){gf[i]=gBus[i].V oltage*sin(gBus[i].Phase);ge[i]=gBus[i].V oltage*cos(gBus[i].Phase);}else{gf[i]=0;ge[i]=1;}}}void GetUnbalance(){int i,j;FILE *fp;for(i=0;i<=Bus_Num-2;i++){gDelta_P[i]=gBus[i+1].GenP-gBus[i+1].LoadP;if(gBus[i+1].Type==2) //PV节点gDelta_Q[i]=gBus[i+1].V oltage*gBus[i+1].V oltage-(ge[i+1]*ge[i+1]+gf[i+1]*gf[i+1]);elsegDelta_Q[i]=gBus[i+1].GenQ-gBus[i+1].LoadQ;for(j=0;j<=Bus_Num-1;j++){gDelta_P[i]=gDelta_P[i]-ge[i+1]*(gY_G[i+1][j]*ge[j]-gY_B[i+1][j]*gf[j])-gf[i+1]*(gY_G[i +1][j]*gf[j]+gY_B[i+1][j]*ge[j]);if(gBus[i+1].Type==1) //PQ节点gDelta_Q[i]=gDelta_Q[i]-gf[i+1]*(gY_G[i+1][j]*ge[j]-gY_B[i+1][j]*gf[j])+ge[i+1]*(gY_G[i+1][ j]*gf[j]+gY_B[i+1][j]*ge[j]);}}for(i=0;i<=Bus_Num-2;i++) //合并{gDelta_PQ[2*i]=gDelta_P[i];gDelta_PQ[2*i+1]=gDelta_Q[i];}if((fp=fopen("F:\\1061180523\\flow\\data\\unbalance.txt","w"))==NULL){printf("Can not open the file named 'unbalance.txt' \n");return ;}fprintf(fp,"---Unbalance---\n");for(i=0;i<=2*Bus_Num-3;i++){fprintf(fp,"Unbalance[%d]=%10.5f\n",i+1,gDelta_PQ[i]);}fclose(fp);}void GetJaccobi(){int i,j;float ga[Bus_Num-1],gb[Bus_Num-1];FILE *fp;for(i=0;i<=Bus_Num-2;i++) //计算注入电流{ga[i]=0;gb[i]=0;for(j=0;j<=Bus_Num-1;j++){ga[i]=ga[i]+gY_G[i+1][j]*ge[j]-gY_B[i+1][j]*gf[j];gb[i]=gb[i]+gY_G[i+1][j]*gf[j]+gY_B[i+1][j]*ge[j];}}for(i=0;i<=Bus_Num-2;i++){for(j=0;j<=Bus_Num-2;j++){if(i!=j){gJaccobi[2*i][2*j]=-gY_B[i+1][j+1]*ge[i+1]+gY_G[i+1][j+1]*gf[i+1];gJaccobi[2*i][2*j+1]=gY_G[i+1][j+1]*ge[i+1]+gY_B[i+1][j+1]*gf[i+1];if(gBus[i+1].Type==2) //PV节点{gJaccobi[2*i+1][2*j]=0;gJaccobi[2*i+1][2*j+1]=0;}else //PQ{gJaccobi[2*i+1][2*j]=-gJaccobi[2*i][2*j+1];gJaccobi[2*i+1][2*j+1]=gJaccobi[2*i][2*j];}}else{gJaccobi[2*i][2*j]=-gY_B[i+1][j+1]*ge[i+1]+gY_G[i+1][j+1]*gf[i+1]+gb[i];gJaccobi[2*i][2*j+1]=gY_G[i+1][j+1]*ge[i+1]+gY_B[i+1][j+1]*gf[i+1]+ga[i];if(gBus[i+1].Type==2) //PV节点{gJaccobi[2*i+1][2*j]=2*gf[i+1];gJaccobi[2*i+1][2*j+1]=2*ge[i+1];}else //PQ{gJaccobi[2*i+1][2*j]=-gY_G[i+1][j+1]*ge[i+1]-gY_B[i+1][j+1]*gf[i+1]+ga[i];gJaccobi[2*i+1][2*j+1]=-gY_B[i+1][j+1]*ge[i+1]+gY_G[i+1][j+1]*gf[i+1]-gb[i];}}}}if((fp=fopen("F:\\1061180523\\flow\\data\\jaccobi.txt","w"))==NULL){printf("Can not open the file named 'jaccobi.txt' \n");return ;}fprintf(fp,"---Jaccobi Matrix---\n");for(i=0;i<=2*Bus_Num-3;i++){for(j=0;j<=2*Bus_Num-3;j++){fprintf(fp,"jaccobi(%d,%d)=%10.5f\n",i+1,j+1,gJaccobi[i][j]);}}fclose(fp);}void GetRevised(){int i,j;FILE *fp;NEquation ob1; //解矩阵方程ob1.SetSize(2*(Bus_Num-1));for(i=0;i<=2*Bus_Num-3;i++)for(j=0;j<=2*Bus_Num-3;j++)ob1.Data(i,j)=gJaccobi[i][j];for(i=0;i<=2*Bus_Num-3;i++)ob1.V alue(i)=gDelta_PQ[i];ob1.Run();for(i=0;i<=Bus_Num-2;i++){gDelta_f[i]=ob1.V alue(2*i);gDelta_e[i]=ob1.V alue(2*i+1);gDelta_fe[2*i]=gDelta_f[i];gDelta_fe[2*i+1]=gDelta_e[i];}if((fp=fopen("F:\\1061180523\\flow\\data\\revised.txt","w"))==NULL) {printf("Can not open the file named 'revised.txt' \n");return ;}fprintf(fp,"---Revised---\n");for(i=0;i<=2*Bus_Num-3;i++){fprintf(fp,"revised[%d]=%10.5f\n",i+1,gDelta_fe[i]);}fclose(fp);}void GetNew V alue(){int i;FILE *fp;for(i=0;i<=Bus_Num-2;i++){gf[i+1]=gf[i+1]+gDelta_f[i];ge[i+1]=ge[i+1]+gDelta_e[i];}if((fp=fopen("F:\\1061180523\\flow\\data\\newvalue.txt","w"))==NULL) {printf("Can not open the file named 'newvalue.txt' \n");return ;}fprintf(fp,"---New V alue---\n");for(i=0;i<=Bus_Num-2;i++){fprintf(fp,"f(%d)=%10.5f,e(%d)=%10.5f\n",i+1,gf[i+1],i+1,ge[i+1]);}fclose(fp);}int main(int argc, char* argv[]){int i,Count_Num;float maxV alue;test();GetData();GetYMatrix();SetInitial();for(Count_Num=0;Count_Num<=100;Count_Num++){GetUnbalance();GetJaccobi();GetRevised();GetNewV alue();maxV alue=fabs(gDelta_fe[0]);for(i=1;i<=2*(Bus_Num-1)-1;i++){if(maxV alue<fabs(gDelta_fe[i])){maxV alue=fabs(gDelta_fe[i]);}}if(maxV alue<Precision){break;}}printf("%d\n",Count_Num);for(i=0;i<=Bus_Num-1;i++){printf("%10.5f\n",sqrt(ge[i]*ge[i]+gf[i]*gf[i]));}return 0;}。

电力系统分析课程设计设计题目9节点电力网络潮流计算指导教师院（系、部）电气与控制工程学院专业班级学号姓名日期电气工程系课程设计标准评分模板目录1 PSASP软件简介 (1)1.1 PSASP平台的主要功能和特点 (1)1.2 PSASP的平台组成 (2)2 牛顿拉夫逊潮流计算简介 (3)2.1 牛顿—拉夫逊法概要 (3)2.2 直角坐标下的牛顿—拉夫逊潮流计算 (5)2.3 牛顿—拉夫逊潮流计算的方法 (6)3 九节点系统单线图及元件数据 (7)3.1 九节点系统单线图 (7)3.2 系统各项元件的数据 (8)4 潮流计算的结果 (10)4.1 潮流计算后的单线图 (10)4.2 潮流计算结果输出表格 (10)5 结论 (14)电力系统分析课程设计任务书9节点系统单线图如下：基本数据如下：母线名基准电压区域号电压上限电压下限发电 1 16.5000 2 18.1500 14.8500 发电 2 18.000 1 19.800 16.2000 发电 3 13.8000 1 15.1800 12.4200 GEN1-230 230.000 2 0.0000 0.0000 GEN2-230 230.000 1 0.0000 0.0000 GEN3-230 230.000 1 0.0000 0.0000 STNA-230 230.000 2 0.0000 0.0000 STNB-230 230.000 2 0.0000 0.0000 STNC-230 230.000 1 0.0000 0.0000数据组I 侧母线J 侧母线编号所属区域单位正序电阻正序电抗正序充电电纳的1/2常规GEN1-230 STNA-230 1 I侧标么0.010000 0.085000 0.088000 常规STNA-230 GEN2-230 2 I侧标么0.032000 0.161000 0.153000 常规GEN2-230 STNC-230 3 I侧标么0.008500 0.072000 0.074500 常规STNC-230 GEN3-230 4 I侧标么0.011900 0.100800 0.104500 常规GEN3-230 STNB-230 5 I侧标么0.039000 0.170000 0.179000 常规STNB-230 GEN1-230 6 I侧标么0.017000 0.092000 0.079000表3 两绕组变压器数据负荷数据1 PSASP软件简介“电力系统分析综合程序”（Power System Analysis Software Package,PSASP）是一套历史悠久、功能强大、使用方便的电力系统分析程序，是高度集成和开发具有我国自主知识产权的大型软件包。

clc;clear;close all;n1=9;n2=9;B1=load('b1.txt','ascii');B2=load('b2.txt','ascii');m=0;r=0;for i=1:n1if B1(i,8)==1m=m+1;endendr=n1-m-1;e=B1(1:n1,2);f=B1(1:n1,3);ps=B1(1:n1,4)-B1(1:n1,6);qs=B1(1:n1,5)-B1(1:n1,7);Y=zeros(n1);for j=1:n2p=B2(j,2);q=B2(j,3);J=zeros(2*n1);Y(p,q)=-1./((B2(j,4)+1i*B2(j,5))*B2(j,7));Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B2(j,4)+1i*B2(j,5))+1i*B2(j,6);Y(q,q)=Y(q,q)+1/((B2(j,4)+1i*B2(j,5))*B2(j,7)^2)+1i*B2(j,6);end%导纳矩阵YH1=zeros(n1);N1=zeros(n1);J1=zeros(n1);L1=zeros(n1);G=real(Y);B=imag(Y);dp(1:(n1-1))=1;dq(1:m)=1;dx(1:(2*m+r))=1;while(abs(max(dp)>=1.0e-5)||abs(max(dq)>=1.0e-5)||abs(max(dx(1:(m+r)))>=1.0e-5)||abs(max(dx((m +r+1):(2*m+r)).*e(1:m))>=1.0e-5))for i=1:n1for j=1:n1si(i,j)=sin(f(i)-f(j));co(i,j)=cos(f(i)-f(j));if i==jH1(i,j)=e(i)^2*B(i,i)+qs(i);N1(i,j)=-e(i)^2*G(i,i)-ps(i);J1(i,j)=e(i)^2*G(i,i)-ps(i);L1(i,j)=e(i)^2*B(i,i)-qs(i);elseH1(i,j)=-e(i)*e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j));N1(i,j)=-e(i)*e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j));J1(i,j)=e(i)*e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j));L1(i,j)=-e(i)*e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j));endendendH=H1(2:n1,2:n1);N=N1(2:n1,r+2:n1);J=J1(r+2:n1,2:n1);L=L1(r+2:n1,r+2:n1);J=[H N;J L];%雅克比矩阵Jfor i=2:n1a(i)=0;b(i)=0;for j=1:n1a(i)=a(i)+e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j));if i>=(r+2)b(i)=b(i)+e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j));endenddp(i-1)=ps(i)-e(i)*a(i);if i>=(r+2)dq(i-(r+1))=qs(i)-e(i)*b(i);endendF=[dp';dq'];dx=-J\F;f(2:n1)=f(2:n1)+dx(1:(m+r));e1=e((r+2):n1)+dx((m+r+1):(2*m+r)).*e((r+2):n1); e2=e(1:(r+1));e=[e2;e1];enddisp('各节点电压');disp('e=');disp(e');disp('单位：度f=');disp(f'*180/3.1415926)%求出各节点电压幅值及相角a=0;for i=1:n1a=a+conj(Y(1,i))*(e(i)*cos(f(i))-1i*e(i)*sin(f(i))); enda=a*(e(1)*cos(f(1))+1i*e(1)*sin(f(1)));disp('平衡节点P+jQ=');disp(a);disp('输电线路有功无功:');for i=1:n2p=B2(i,2);q=B2(i,3);s1(p,q)=e(p)^2*(-1i*B2(i,6))+(e(p)*cos(f(p))+1i*e(p)*sin(f(p)))*(e(p)*cos(f(p))-1i*e(p)*sin(f(p))-e( q)*cos(f(q))+1i*e(q)*sin(f(q)))*conj(-Y(q,p));s1(q,p)=e(q)^2*(-1i*B2(i,6))+(e(q)*cos(f(q))+1i*e(q)*sin(f(q)))*(e(q)*cos(f(q))-1i*e(q)*sin(f(q))-e( p)*cos(f(p))+1i*e(p)*sin(f(p)))*conj(-Y(q,p));fprintf('S%d%d=',p,q);disp(s1(p,q));fprintf('S%d%d=',q,p);disp(s1(q,p));end%初值计算for j=1:n1s(j)=0;for i=1:n1s(j)=s(j)+conj(Y(j,i))*(e(i)*cos(f(i))-1i*e(i)*sin(f(i)));ends(j)=s(j)*(e(j)*cos(f(j))+1i*e(j)*sin(f(j)));endY1=conj(-s)./(e'.^2);%负荷等值导纳disp('各节点等值导纳');disp('Y1=');disp(Y1);disp('各节点功率');disp('S=');disp(s);Ra=[0 0 0];X_d=[0.0608 0.1198 0.1813];Xq=[0.0969 0.8645 1.2578];%Xq=X_d;%不计凸极效应v=e.*cos(f)+1i*e.*sin(f);for i=1:3I(i)=conj(s(i)/(e(i)*cos(f(i))+1i*e(i)*sin(f(i))));EQ(i)=e(i)*cos(f(i))+1i*e(i)*sin(f(i))+(Ra(i)+1i*Xq(i))*I(i);EQx(i)=real(EQ(i));EQy(i)=imag(EQ(i));delta(i)=atan(EQy(i)/EQx(i));enddisp('各发电机的暂态电动势，功角和输入机械功率初值');disp('delta=');disp(delta*180/3.1415926);for i=1:3Vx(i)=e(i)*cos(f(i));Vy(i)=e(i)*sin(f(i));Vdq=[sin(delta(i)) -cos(delta(i));cos(delta(i)) sin(delta(i))]*[Vx(i);Vy(i)];Vd(i)=Vdq(1);Vq(i)=Vdq(2);Ix(i)=real(I(i));Iy(i)=imag(I(i));Idq=[sin(delta(i)) -cos(delta(i));cos(delta(i)) sin(delta(i))]*[Ix(i);Iy(i)];Id(i)=Idq(1);Iq(i)=Idq(2);E_q(i)=Vq(i)+Ra(i)*Iq(i)+X_d(i)*Id(i);enddisp('E_q=');disp(E_q);for i=1:3Pm(i)=real(s(i))+(Ix(i)^2+Iy(i)^2)*Ra(i); enddisp('Pm=');disp(Pm);。

#include <iostream>#include <fstream>#include<iomanip>#include<math.h>using namespace std;//节点号类型负荷有功负荷无功母线数据（类型1＝PV节点，2＝PQ节点，3＝平衡节点）struct BUS{int busno;int type;float Pd;float Qd;};//发电机数据节点号有功发电电压幅值struct Generator{int busno;float Pg;float Vg;};//支路信息节点I 节点J R X B/2 kstruct Line{int busi;int busj;float R;float X;float B;float k;};//deltaP deltaQ deltaV^2//void fun1(double YG[][50],double YB[][50],double e[],double f[],int type[],int N,double W[],double P[],double Q[],double V[]){double dP=0,dQ=0,dV=0;int i,j;for(i=0;i<N-1;i++){double A=0,B=0;for(j=0;j<N;j++){A+=YG[i][j]*e[j]-YB[i][j]*f[j];B+=YG[i][j]*f[j]+YB[i][j]*e[j];}dV=V[i]*V[i]-e[i]*e[i]-f[i]*f[i];dP=P[i]-e[i]*A-f[i]*B;W[2*i]=dP;dQ=Q[i]-f[i]*A+e[i]*B;if(type[i]==1)W[2*i+1]=dQ;else W[2*i+1]=dV;}}//Jacobi矩阵//void Jacobi(double YG[][50],double YB[][50],double e[50],double f[50],int type[50],int N ,double Ja[100][101]){ int i,j;for(i=0;i<N;i++){for(j=0;j<N;j++){if(i!=j){if(type[i]==1){Ja[2*i][2*j]=-(YG[i][j]*e[i]+YB[i][j]*f[i]);Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i];Ja[2*i+1][2*j]=Ja[2*i][2*j+1];Ja[2*i+1][2*j+1]=-Ja[2*i][2*j];}else {Ja[2*i][2*j]=-YG[i][j]*e[i]+YB[i][j]*f[i];Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i];Ja[2*i+1][2*j+1]=Ja[2*i+1][2*j]=0;}}else {double a[50]={0},b[50]={0};for(int k=0;k<N;k++){a[i]+=(YG[i][k]*e[k]-YB[i][k]*f[k]);b[i]+=(YG[i][k]*f[k]+YB[i][k]*e[k]);Ja[2*i][2*j]=-a[i]-YG[i][i]*e[i]-YB[i][i]*f[i];Ja[2*i][2*j+1]=-b[i]+YB[i][i]*e[i]-YG[i][i]*f[i];if(type[i]==1){Ja[2*i+1][2*j]=b[i]+YB[i][i]*e[i]-YG[i][i]*f[i];Ja[2*i+1][2*j+1]=-a[i]+YG[i][i]*e[i]+YB[i][i]*f[i];}else {Ja[2*i+1][2*j]=-2*e[i];Ja[2*i+1][2*j+1]=-2*f[i];}}}}}}//高斯消元法解方程组函数//void gauss(double a[][101],int n){int i,j,k;double c;for(k=0;k<n-1;k++) {c=a[k][k];for(j=k;j<=n;j++) a[k][j]/=c;for(i=k+1;i<n;i++) {c=a[i][k];for(j=k;j<=n;j++) a[i][j]-=c*a[k][j];}}a[n-1][n]/=a[n-1][n-1];for(k=n-2;k>=0;k--)for(j=k+1;j<n;j++) a[k][n]-=a[k][j]*a[j][n];}void main(){ifstream fin;int N=0,GS=0,LD=0,ZLs=0; //节点数发电机数负荷数支路数// BUS *B;Generator *G;Line *L;//从文本中读入原始数据到数组中//fin.open("C:\\data.txt");if(!fin){cout<<"输入数据文件不存在！"<<endl;getchar();}int m1[50]={0},m2[50]={0};float m3[50],m4[50],m5[50],m6[50];int i,j,l;for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m2[i]>>m3[i]>>m4[i];N++;}B =new BUS[N];for (i=0;i<N;i++){B[i].busno=m1[i];B[i].type=m2[i];B[i].Pd=m3[i];B[i].Qd=m4[i];}for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m4[i]>>m3[i];GS++;}G =new Generator[GS];for (i=0;i<GS;i++){G[i].busno=m1[i];G[i].Pg=m4[i];G[i].Vg=m3[i];}for(i=0;;i++){fin>>m1[i];if(m1[i]==0)break;fin>>m2[i]>>m3[i]>>m4[i]>>m5[i]>>m6[i];ZLs++;}L =new Line[ZLs];for (i=0;i<ZLs;i++){L[i].busi=m1[i];L[i].busj=m2[i];L[i].R=m3[i];L[i].X=m4[i];L[i].B=m5[i];L[i].k=m6[i];}LD=N-GS;fin.close();//节点导纳矩阵形成//double YB[50][50],YG[50][50],BB[50][50],K[50][50];for(i=0;i<N;i++){for(j=0;j<N;j++){YB[i][j]=0;YG[i][j]=0;BB[i][j]=0;K[i][j]=1;}}for (l=0;l<ZLs;l++){i=L[l].busi-1;j=L[l].busj-1;K[i][j]=L[l].k;BB[i][j]=BB[j][i]=L[l].B;YG[i][j]=YG[j][i]=L[l].R/(L[l].R*L[l].R+L[l].X*L[l].X);YB[i][j]=YB[j][i]=-L[l].X/(L[l].R*L[l].R+L[l].X*L[l].X);}for(i=0;i<N;i++){for(j=i;j<N;j++) {K[i][j]=K[j][i];K[j][i]=1;}for(j=0;j<N;j++){if(i!=j){YG[i][i]=YG[i][i]+(YG[i][j]*K[i][j]*K[i][j]);YB[i][i]=YB[i][i]+(YB[i][j]*K[i][j]*K[i][j]+BB[i][j]);}}}//修正后//for (l=0;l<ZLs;l++){i=L[l].busi-1;j=L[l].busj-1;K[i][j]=L[l].k;YG[i][j]=-YG[i][j]*K[i][j];YG[j][i]=YG[i][j];YB[i][j]=-YB[i][j]*K[i][j];YB[j][i]=YB[i][j];}int type[50]={0};for(i=0;i<N;i++){type[i]=B[i].type;}//PQV的获得//double P[50],Q[50],V[50];for(i=0;i<N;i++){P[i]=0;Q[i]=0;V[i]=0;P[i]=-B[i].Pd;Q[i]=-B[i].Qd;}for (i=0;i<GS;i++){P[G[i].busno-1]=G[i].Pg;V[G[i].busno-1]=G[i].Vg;}// 求A=e+f//double e[50]={0},f[50]={0};double C[100]={0},D[100]={0};for(i=0;i<N;i++){if(V[i]==0){C[2*i]=1;}else C[2*i]=V[i];}double W[100]={0},Ja[100][101]={0}; //调用Jacobi函数和高斯函数//for(int t=1;t<10;t++){for(i=0;i<2*N-2;i++){e[i]=C[2*i];f[i]=C[2*i+1];}fun1(YG,YB,e,f,type,N,W,P,Q,V);double it=fabs(W[0]);for(i=1;i<2*N-2;i++){if (it<fabs(W[i])) {it=fabs(W[i]);j=i;}}//中间迭代过程//cout<<setw(10)<<"迭代次数"<<setw(20)<<"最大的功率误差"<<setw(8)<<"节点号"<<endl;cout<<setw(10)<<t<<setw(20)<<it<<setw(8)<<j/2+1<<endl;if (it<0.00001) break;Jacobi(YG,YB,e,f,type,N,Ja);for(i=0;i<2*N-2;i++){Ja[i][2*N-2]=W[i];}//高斯消元法解方程//gauss(Ja,2*N-2);for(i=0;i<2*N-2;i++){D[i]=-Ja[i][2*(N-1)];C[i]+=D[i];}}//平衡节点//for(i=0;i<N;i++){double a=0,b=0;for(int j=0;j<N;j++){a+=(YG[i][j]*e[j]-YB[i][j]*f[j]);b+=(YB[i][j]*e[j]+YG[i][j]*f[j]);}P[i]=e[i]*a+f[i]*b;Q[i]=f[i]*a-e[i]*b;}//支路//double PZL[100][101]={0},QZL[100][101]={0},pr[100][101]={0},qx[100][101]={0};double x1=0,x2=0,y1=0,y2=0,I2=0;for(int k=0;k<ZLs;k++){i=L[k].busi-1;j=L[k].busj-1;x1=e[i]/L[k].k-e[j];y1=f[i]/L[k].k-f[j];x2=-e[i]*YG[i][j]-f[i]*YB[i][j];y2=-f[i]*YG[i][j]+e[i]*YB[i][j];QZL[i][j]=(x1*y2-x2*y1);PZL[i][j]=(x1*x2+y1*y2);I2=(PZL[i][j]*PZL[i][j]+QZL[i][j]*QZL[i][j])/(e[i]*e[i]+f[i]*f[i]);pr[i][j]=I2*L[k].R;qx[i][j]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[i][j]+=(e[i]*e[i]+f[i]*f[i])*(-L[k].B);x1=e[j]*L[k].k-e[i];y1=f[j]*L[k].k-f[i];x2=-e[j]*YG[j][i]-f[j]*YB[j][i];y2=-f[j]*YG[j][i]+e[j]*YB[j][i];QZL[j][i]=(x1*y2-x2*y1);PZL[j][i]=(x1*x2+y1*y2);I2=(PZL[j][i]*PZL[j][i]+QZL[j][i]*QZL[j][i])/(e[j]*e[j]+f[j]*f[j]);pr[j][i]=I2*L[k].R;qx[j][i]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[j][i]+=(e[j]*e[j]+f[j]*f[j])*(-L[k].B);}//全网数据//int high=1,low=1;double PG=0,PL=0,Prr=0,Vh=sqrt(e[0]*e[0]+f[0]*f[0]),Vl=sqrt(e[0]*e[0]+f[0]*f[0]);for(k=0;k<N;k++){Prr+=P[k];if(B[k].type==1) PL+=B[k].Pd;else PG+=P[k];if(sqrt(e[k]*e[k]+f[k]*f[k])>Vh){Vh=sqrt(e[k]*e[k]+f[k]*f[k]);high=k+1;}if(sqrt(e[k]*e[k]+f[k]*f[k])<Vl){Vl=sqrt(e[k]*e[k]+f[k]*f[k]);low=k+1;}}//输出数据到文件databak.txt//ofstream fout;fout.open("C:\\databak.txt");fout<<"节点"<<endl;fout<<setw(8)<<"节点号"<<setw(16)<<"V"<<setw(16)<<"弧度"<<setw(16)<<"发电P"<<setw(16)<<"发电Q"<<setw(16)<<"负荷P"<<setw(16)<<"负荷Q"<<endl;for(i=0;i<LD;i++){fout<<setw(8)<<i+1<<setw(16)<<sqrt(e[i]*e[i]+f[i]*f[i])<<setw(16)<<atan2(f[i],e[i])*180/3.14159<<setw(16)<<0 <<setw(16)<<0<<setw(16)<<B[i].Pd<<setw(16)<<B[i].Qd<<endl;}for(j=0;j<GS;j++){i=G[j].busno-1;fout<<setw(8)<<i+1<<setw(16)<<V[i]<<setw(16)<<atan2(f[i],e[i])*180/3.14159<<setw(16)<<P[i]<<setw(16)<<Q[i]<< setw(16)<<0<<setw(16)<<0<<endl;}fout<<"支路"<<endl;fout<<setw(4)<<"i"<<setw(4)<<"j"<<setw(10)<<"i_j有功"<<setw(10)<<"i_j无功"<<setw(10)<<"j_i有功"<<setw(12)<<"j_i无功"<<setw(12)<<"有功损耗"<<setw(12)<<"无功损耗"<<endl;for (k=0;k<ZLs;k++){i=L[k].busi-1;j=L[k].busj-1;fout<<setw(4)<<L[k].busi<<setw(4)<< L[k].busj <<setw(10)<<PZL[i][j] <<setw(10)<<QZL[i][j]<<setw(10)<< PZL[j][i] <<setw(12)<<QZL[j][i]<<setw(12)<<pr[i][j]<<setw(12)<<qx[i][j]<<endl;}fout<<"全网数据"<<endl;fout<<setw(14)<<"发电有功"<<setw(14)<<"负荷有功"<<setw(14)<<"有功损耗"<<setw(14)<<"最高电压"<<setw(14)<<"节点号"<<setw(14)<<"最低电压"<<setw(14)<<"节点号"<<endl;fout<<setw(14)<<PG<<setw(14)<<PL<<setw(14)<<Prr<<setw(14)<<Vh<<setw(14)<<high<<setw(14)<<Vl<<setw(14) <<low<<endl;fout.close();}。