电力系统潮流计算C语言程序及说明知识分享

电力系统潮流计算的计算机程序设计一、本文概述随着电力系统的日益复杂化和智能化，对电力系统潮流计算的精度和效率提出了更高的要求。

潮流计算作为电力系统分析的基础，其准确性直接关系到电力系统的安全、稳定和经济运行。

本文旨在探讨电力系统潮流计算的计算机程序设计，以提高计算效率，优化计算结果，为电力系统的规划、设计、运行和管理提供有力支持。

本文首先介绍了电力系统潮流计算的基本原理和方法，包括节点导纳矩阵的形成、功率方程的求解等。

在此基础上，详细阐述了潮流计算的计算机程序设计，包括程序设计的总体思路、主要模块的功能和实现方法。

同时，结合具体的算例和仿真实验，对程序设计的有效性进行了验证和分析。

本文还讨论了潮流计算程序设计中的关键技术和难点，如数值稳定性、收敛性等问题，并提出了相应的解决策略。

还对潮流计算程序设计的未来发展趋势进行了展望，包括考虑更多约束条件、引入智能优化算法、实现并行计算等方面的研究和应用。

本文旨在通过计算机程序设计的角度，深入探讨电力系统潮流计算的理论和实践，为电力系统的安全运行和可持续发展提供有益的技术支持和指导。

二、电力系统基础知识电力系统是指由发电、输电、变电、配电和用电等环节组成的电能生产与消费系统。

它不仅是保证电能生产、输送、分配和使用的系统，也是一个庞大而复杂的工程系统。

在电力系统中，潮流计算是一项至关重要的任务，它决定了电网的运行状态，为电力系统的稳定、经济、安全运行提供了重要依据。

电力系统的基本构成主要包括发电厂、输电线路、变压器、配电线路和用户。

发电厂负责将一次能源转化为电能，输电线路负责将电能从发电厂输送到各个变电站，变压器则负责调整电压等级以满足不同用户的需求，配电线路则将电能从变电站输送到各个用户，而用户则是电能的最终消费者。

在电力系统中，电压和电流是描述电能状态的两个基本物理量。

电压是指电场中单位正电荷移动的势能差，通常用字母U或V表示。

电流则是指单位时间内通过导体横截面的电荷量，通常用字母I表示。

程序的稳定性分析
程序在不同计算机上的运行 结果是否一致。
程序运行过程中，输入数据 的变化对输出结果的影响程 度。
程序在长时间运行过程中， 输出结果是否保持稳定。
程序在处理异常情况时的表 现和稳定性。
程序的扩展性分析
代码可读性：C语言程序应具备良好的可读性，方便后续维护和修改 算法效率：C语言程序应采用高效的算法，提高计算速度 内存占用：C语言程序应合理利用内存，避免内存泄漏和不必要的内存占用 扩展性：C语言程序应具备良好的扩展性，方便添加新功能和优化性能
THANK YOU
汇报人：XX
程序的异常处理说明
异常类型：输入 错误、计算错误、 内存不足等
异常处理方式： 使用try-catch 语句进行异常捕 获和处理
异常处理流程： 当异常发生时， 程序会输出错误 信息并终止运行
异常处理结果： 确保程序在遇到 异常时能够正确 处理并给出相应 的提示信息
C语言程序应用示例
示例程序的输入数据格式
添加标题
添加标题
添加标题Βιβλιοθήκη 输入输出函数：用于数据的输入和 输出
函数：可重复使用的代码块，具有 特定的功能
C语言程序中电力系统模型的建立
定义节点和支路：根 据电力系统网络结构， 定义节点和支路，为 潮流计算做准备。
建立数学模型：根据 电力系统的物理特性 和元件参数，建立数 学模型，包括节点电 压、支路电流等。
实际运行时 间测试
程序的内存占用性能分析
内存占用情况：分 析程序运行过程中 内存的占用情况， 包括堆内存和栈内 存。
内存泄漏检测：检 查程序是否存在内 存泄漏，即程序运 行结束后未正确释 放的内存。
内存使用优化：根 据内存占用情况， 优化程序中的数据 结构或算法，降低 内存占用。

for(j=0;j<Bus_Num;j++)
{
ia[i]=ia[i]+gY_G[n][j]*ge[j]-gY_B[n][j]*gf[j];
ib[i]=ib[i]+gY_G[n][j]*gf[j]+gY_B[n][j]*ge[j];
}
}
for(i=0,n=1;i<Bus_Num-1;i++,n++)
{
gDelta_PQ[2*i]=gDelta_P[i];
gDelta_PQ[2*i+1]=gDelta_Q[i];
}
if((fp=fopen("C:\\Documents and Settings\\Zorro\\桌面\\1\\data\\unbalance.txt","w"))==NULL)
if(gBus[n].Type==1)
gDelta_Q[i]=gDelta_Q[i]-gf[n]*(gY_G[n][j]*ge[j]-gY_B[n][j]*gf[j])+ge[n]*(gY_G[n][j]*gf[j]+gY_B[n][j]*ge[j]);
}
}
for(i=0;i<Bus_Num-1;i++)
{
gY_G[i][j]=0.0;
gY_B[i][j]=0.0;
}
for(l=0;l<Line_Num;l++)
{
i=gLine[l].No_I-1;
j=gLine[l].No_J-1;
r=gLine[l].R;
x=gLine[l].X;

班级：姓名：学号：一、作业要求编写程序计算图1所示算例系统的潮流及三相短路电流..潮流计算：方法不限;计算系统的节点电压和相角..短路电流：4号母线发生金属性三相短路时z f=0;分别按照精确算法和近似算法计算短路电流、系统中各节点电压以及网络中各支路的电流分布;并对两种情况下的计算结果进行比较..二、电路图及参数图1 3机9节点系统表1 9节点系统支路参数表2 9节点系统发电机参数表3 9节点系统负荷参数三、计算步骤(1) 进行系统正常运行状态的潮流计算;求得(0)i U (2) 形成不含发电机和负荷的节点导纳矩阵Y N ；(3) 将发电机表示为电流源i I /i diE jx ''=和导纳i y 1/di jx '=的并联组合；节点负荷用恒阻抗的接地支路表示；形成包括所有发电机支路和负荷支路的节点导纳矩阵Y;即在Y N 中的发电机节点和负荷节点的自导纳上分别增加发电机导纳i y 和负荷导纳,LD i y *,,22LD i LDi LDiLD i i i S P jQ y V V -==； (4) 利用1Z Y-=;计算节点阻抗矩阵;从而得到阻抗矩阵中的第f 列；(5) 利用公式6-7或6-10计算短路电流；(6)利用公式6-8或6-11计算系统中各节点电压；(7)利用公式6-9计算变压器支路的电流；对输电线路利用П型等值电路计算支路电流..四、计算结果节点导纳矩阵Yn：Columns 1 through 50 -17.3611i 0 0 0 +17.3611i 00 0 -16.0000i 0 0 00 0 0 -17.0648i 0 00 +17.3611i 0 0 3.3074 -39.3089i -1.3652 +11.6041i0 0 0 -1.3652 +11.6041i 2.5528 -17.3382i0 0 0 -1.9422 +10.5107i 00 0 +16.0000i 0 0 -1.1876 + 5.9751i0 0 0 0 00 0 0 +17.0648i 0 0 Columns 6 through 90 0 0 00 0 +16.0000i 0 00 0 0 0 +17.0648i-1.9422 +10.5107i 0 0 00 -1.1876 + 5.9751i 0 03.2242 -15.8409i 0 0 -1.2820 + 5.5882i0 2.8047 -35.4456i -1.6171 +13.6980i 00 -1.6171 +13.6980i 2.7722 -23.3032i -1.1551 + 9.7843i-1.2820 + 5.5882i 0 -1.1551 + 9.7843i 2.4371 -32.1539i电压幅值：1.0400 1.0250 1.0250 1.0258 0.9956 1.0127 1.0258 1.0159 1.0324电压相角：0 0.1620 0.0814 -0.0387 -0.0696 -0.0644 0.0649 0.0127 0.0343节点有功：0.7164 1.6300 0.8500 0.0000 -1.2500 -0.9000 -0.0000 -1.0000 -0.0000节点无功：0.2705 0.0665 -0.1086 0.0000 -0.5000 -0.3000 -0.0000 -0.3500 -0.0000修正后的节点导纳矩阵Y：Columns 1 through 50 -20.6944i 0 0 0 +17.3611i 00 0 -19.3333i 0 0 00 0 0 -20.3982i 0 00 +17.3611i 0 0 3.3074 -39.3089i -1.3652 +11.6041i0 0 0 -1.3652 +11.6041i 3.8716 -17.6627i0 0 0 -1.9422 +10.5107i 00 0 +16.0000i 0 0 -1.1876 + 5.9751i0 0 0 0 00 0 0 +17.0648i 0 0 Columns 6 through 90 0 0 00 0 +16.0000i 0 00 0 0 0 +17.0648i-1.9422 +10.5107i 0 0 00 -1.1876 + 5.9751i 0 04.1321 -16.0184i 0 0 -1.2820 +5.5882i0 2.8047 -35.4456i -1.6171 +13.6980i 00 -1.6171 +13.6980i 3.7323 -23.6669i -1.1551 + 9.7843i-1.2820 + 5.5882i 0 -1.1551 + 9.7843i 2.4371 -32.1539i节点阻抗矩阵Z的第4列：0.0463 + 0.1252i0.0329 + 0.0693i0.0316 + 0.0707i0.0552 + 0.1493i0.0589 + 0.1204i0.0562 + 0.1226i0.0397 + 0.0838i0.0416 + 0.0814i0.0378 + 0.0845i精确计算结果：短路电流：模值：6.4459相角：-71.9365节点电压模值：0.1831 0.5687 0.5427 0.0000 0.1466 0.1506 0.4537 0.4463 0.4495支路电流：i j Iij1 4 0.5779-3.1264i2 7 1.3702-1.4433i3 9 0.64294-1.4808i4 5 -0.77968+1.5248i4 6 -0.6411+1.477i5 7 -0.89528+1.6436i6 9 -0.73353+1.5487i7 8 0.50734+0.10234i8 9 0.062766+0.056451i近似计算结果：短路电流：模值：6.2838相角：-69.7198节点电压模值：0.1611 0.5214 0.5157 0.0000 0.1827 0.1675 0.4227 0.4348 0.4217五、程序流程图六、程序及输入文件input_data.xls 文件：powerflow_cal.m 文件：l=9;%支路数n=9;%节点数m=6;%PQ节点数Yn=zerosn;%初始化节点导纳矩阵Y DATA1=xlsread'input_data.xls';1; %计算节点导纳矩阵Yfor k=1:li=DATA1k;1;j=DATA1k;2;R=DATA1k;3;X=DATA1k;4;B2=DATA1k;5;Yni;i=Yni;i+1i*B2+1/R+1i*X; Ynj;j=Ynj;j+1i*B2+1/R+1i*X; Yni;j=Yni;j-1/R+1i*X;Ynj;i=Ynj;i-1/R+1i*X;enddisp'节点导纳矩阵Yn：';dispYn;G=realYn;B=imagYn;DATA2=xlsread'input_data.xls';2;P=zeros1;n;Q=zeros1;n;U=ones1;n;P2:n=DATA22:n;3;Q4:n=DATA24:n;4;U1:3=DATA21:3;5;%设置节点电压初值e1=DATA21;5;e2:n=1.0;f1:n=0.0;%设置迭代次数t=0;tmax=10;while t<=tmax%计算fxa1:n=0.0;c1:n=0.0;for i=2:nfor j=1:nai=ai+Gi;j*ej-Bi;j*fj;ci=ci+Gi;j*fj+Bi;j*ej;endendfor i=2:ndeltaPi=Pi-ei*ai-fi*ci;endfor j=4:ndeltaQj=Qj-fj*aj+ej*cj;endfor k=2:3deltaU2k=Uk*Uk-ek*ek-fk*fk;endfx=deltaP2:n deltaQ4:n deltaU22:3';%计算雅克比矩阵Jfor i=2:nfor j=2:nif i~=jHi;j=-Gi;j*ei+Bi;j*fi;Ni;j=Bi;j*ei-Gi;j*fi;elseHi;j=-ai-Gi;i*ei+Bi;i*fi;Ni;j=-ci+Bi;i*ei-Gi;i*fi;endendendfor i=4:nfor j=2:nif i~=jMi;j=Bi;j*ei-Gi;j*fi;Li;j=Gi;j*ei+Bi;j*fi;elseMi;j=ci+Bi;i*ei-Gi;i*fi;Li;j=-ai+Gi;i*ei+Bi;i*fi;endendendfor i=2:3for j=2:nif i~=jRi;j=0;Si;j=0;elseRi;j=-2*ei;Si;j=-2*fi;endendendJ=H2:n;2:n N2:n;2:n;M4:n;2:n L4:n;2:n;R2:3;2:n S2:3;2:n;if maxabsfx<0.0001%输出结果break;else%求解修正方程获得dxdx=-J^-1*fx;dx=dx';e2:n=e2:n+dx1:n-1;f2:n=f2:n+dxn:2*n-1;t=t+1;endendif t>tmaxstr='潮流计算不收敛';dispstr;elsea1:n=0.0;c1:n=0.0;for i=1:nfor j=1:nai=ai+Gi;j*ej-Bi;j*fj; ci=ci+Gi;j*fj+Bi;j*ej;endendfor i=1:nUi=ei+1i*fi;ampi=absUi;argi=angleUi;Pi=ei*ai+fi*ci;Qi=fi*ai-ei*ci;enddisp'电压幅值：';dispamp;disp'电压相角：';disparg;disp'节点有功：';dispP;disp'节点无功：';dispQ;end%计算短路电流f=4;zf=0.0;%修正节点导纳矩阵Xd=DATA21:3;6;E=DATA21:3;7;for i=1:3Iii=Ei/1i*Xdi;endY=Yn;for i=1:3Yi;i=Yi;i+1/1i*Xdi;endfor j=4:nYj;j=Yj;j+-Pj+1i*Qj/Uj*Uj;enddisp'修正后的节点导纳矩阵Y：';Z=Y^-1;disp'节点阻抗矩阵Z的第4列：';dispZ:;4;%精确计算disp'精确计算结果：';U0=U;If=U0f/Zf;f+zf;amp=absIf;arg=atandimagIf/realIf;disp'短路电流：';disp'模值：';dispamp;disp'相角：';disparg;for i=1:nUi=U0i-Zi;f*If;amp=absU;enddisp'节点电压模值：';dispamp;disp'支路电流： ';str='i ''j '' Iij';dispstr;for k=1:li=DATA1k;1;j=DATA1k;2;r=DATA1k;3;x=DATA1k;4;z=r+1i*x;I=Ui-Uj/z;str=num2stri ' ' num2strj ' ' num2strI; dispstr;end%近似计算disp'近似计算结果：';U01:n=1.0;If=U0f/Zf;f+zf;amp=absIf;arg=atandimagIf/realIf;disp'短路电流：';disp'模值：';dispamp;disp'相角：';for i=1:nUi=U0i-Zi;f*If; amp=absU;enddisp'节点电压模值：'; dispamp;。

电力系统潮流计算完整程序及详细理论说明——秦羽风在我刚开始学习潮流程序时,总是找不到一个正确的程序开始模仿学习。

后来经过多方努力，终于自己写出了一个结构清晰、完整的潮流程序。

此程序是一个通用的程序，只需要修改输入数据的子函数（PowerFlowsData_K)里面的母线、支路、发电机、负荷，就能算任意一个网络结构的交流系统潮流。

很适合初学者学习.为了帮助电力系统的同学一起学习,我将我编写的潮流计算程序分享下来给大家；此程序是在基于牛顿拉夫逊算法的基础上，编写的快速解耦算法。

每一个子程序我都有备注说明。

如果有不对的地方，希望大家指正！下文中呈现的顺序为：网络结构、子程序、主程序、运算结果、程序设计理论说明。

一、网络结构：5节点网络如下图。

二、子程序（共有9个子程序）子程序1：（其他系统，只需要修改Bus、Branch、Generator、Load,这四个矩阵就行了）function [Bus,Branch，Generator,Load］=PowerFlowsData_K％％节点数据% 类型:1-平衡节点；2-发电机PV节点；3—负荷PQ节点；4-发电机PQ节点;Bus=［% 类型电压相角1 1。

06 0；2 1 0；3 1 0;3 1 0；3 1 0]；％% 线路数据Branch=[% 发送接收电阻电感(电导电容)并联1 2 0.02 0.06 0 0.06；1 3 0。

08 0。

24 0 0。

05；2 3 0.06 0.18 0 0。

04;2 4 0。

06 0。

18 0 0.04;2 5 0.04 0.12 0 0。

03；3 4 0.01 0.03 0 0。

02;4 5 0.08 0.24 0 0.05］;%% 发电机数据Generator=［％节点定有功定无功(上限下限）无功1 0 0 5 —5;2 0。

4 03 —3];％％负载数据Load=［% 节点定有功定无功2 0.2 0.1；3 0。

【最新整理，下载后即可编辑】电力系统潮流计算注：这是一个基于N－R法的潮流计算通用程序，仅提供了子程序，需要做些处理才能成为一个可运行的计算程序！此程序非我原创，仅与大家共享！！！/*************************************************************** ***** 这里提供的是电力系统潮流计算机解法的五个子程序,采用的方法是** Newton_Raphson法. ** 程序中所用的变量说明如下: ** N:网络节点总数. M:网络的PQ节点数. ** L:网络的支路总数. N0:雅可比矩阵的行数. ** N1:N0+1 K:打印开关.K=1,则打印;否则,不打印.** K1:子程序PLSC中判断输入电压的形式.K1=1,则为极座标形式.否则** 为直角坐标形式. ** D:有功及无功功率误差的最大值. * * G(I,J):Ybus的电导元素(实部). ** B(I,J):Ybus的电纳元素(虚部). ** G1(I) :第I支路的串联电导. B1(I):第I支路的串联电纳. ** C1(I) :第I支路的pie型对称接地电纳. * * C(I,J):第I节点J支路不对称接地电纳. * * CO(I) :第I节点的接地电纳. ** S1(I) :第I节点的起始节点号. E1(I):第I节点的终止节点号. ** P(I) :第I节点的注入有功功率. Q(I):第I节点的注入无功功率.** P0(I) :第I节点有功功率误差. Q0(I):第I节点无功功率误差. ** V0(I) :第I节点(PV节点)的电压误差(平方误差). ** V(I) :第I节点的电压误差幅值. * * E(I) :第I节点的电压的实部. F(I):第I节点的电压的虚部. ** JM(I,J):Jacoby矩阵的第I行J列元素. * * A(I,J):修正方程的增广矩阵,三角化矩阵的第I行J列元素,运算结** 束后A矩阵的最后一列存放修正的解. ** P1(I) :第I支路由S1(I)节点注入的有功功率. ** Q1(I) :第I支路由S1(I)节点注入的无功功率. ** P2(I) :第I支路由E1(I)节点注入的有功功率. ** Q2(I) :第I支路由E1(I)节点注入的无功功率. ** P3(I) :第I支路的有功功率损耗. ** Q3(I) :第I支路的无功功率损耗. ** ANGLE(I):第I节点电压的角度. * **************************************************************** ***/#include <math.h>#include <stdio.h>#define f1(i) (i-1)/* 把习惯的一阶矩阵的下标转化为C语言数组下标*/#define f2(i,j,n) ((i-1)*(n)+j-1)/* 把习惯的二阶矩阵的下标转化为C语言数组下标*//***************************************************** 本子程序根据所给的支路导纳及有关信息,形成结点** 导纳矩阵,如打印参数K=1,则输出电导矩阵G和电纳矩B * ****************************************************/void ybus(int n,int l,int m,float *g,float *b,float *g1,float *b1,float *c1,\ float *c,float *co,int k,int *s1,int *e1){extern FILE *file4;FILE *fp;int i,j,io,i0;int pos1,pos2;int st,en;if(file4==NULL){fp=stdout;}else{fp=file4; /* 输出到文件*/}/* 初始化矩阵G,B */for(i=1;i<=n;i++){for(j=1;j<=n;j++){pos2=f2(i,j,n);g[pos2]=0;b[pos2]=0;}}/* 计算支路导纳*/for(i=1;i<=l;i++){/* 计算对角元*/pos1=f1(i);st=s1[pos1];en=e1[pos1];pos2=f2(st,st,n);g[pos2]+=g1[pos1];b[pos2]+=b1[pos1]+c1[pos1];pos2=f2(en,en,n);g[pos2]+=g1[pos1];b[pos2]+=b1[pos1]+c1[pos1];/* 计算非对角元*/pos2=f2(st,en,n);g[pos2]-=g1[pos1];b[pos2]-=b1[pos1];g[f2(en,st,n)]=g[f2(st,en,n)];b[f2(en,st,n)]=b[f2(st,en,n)]; }/* 计算接地支路导纳*/for(i=1;i<=n;i++){/* 对称部分*/b[f2(i,i,n)]+=co[f1(i)];/* 非对称部分*/for(j=1;j<=l;j++){b[f2(i,i,n)]+=c[f2(i,j,l)];}}if(k!=1){return; /* 如果K不为1,则返回;否则,打印导纳矩阵*/}fprintf(fp,"\n BUS ADMITTANCE MATRIX Y(BUS):"); fprintf(fp,"\n ******************* ARRAY G ********************"); for(io=1;io<=n;io+=5){i0=(io+4)>n?n:(io+4);fprintf(fp,"\n");for(j=io;j<=i0;j++){fprintf(fp,"%13d",j);}for(i=1;i<=n;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%13.6f",g[f2(i,j,n)]);}}fprintf(fp,"\n");}fprintf(fp,"\n ******************* ARRAY B ********************"); for(io=1;io<=n;io+=5){i0=(io+4)>n?n:(io+4);fprintf(fp,"\n");for(j=io;j<=i0;j++){fprintf(fp,"%13d",j);}for(i=1;i<=n;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%13.6f",b[f2(i,j,n)]);}}fprintf(fp,"\n");}fprintf(fp,"\n************************************************"); }/******************************************** 本子程序根据所给的功率及电压等数据** 求出功率及电压误差量,并返回最大有功功率** 以用于与给定误差比较.如打印参数K=1,则输** 出P0,Q0(对PQ结点),V0(对PV结点). ********************************************/void dpqc(float *p,float *q,float *p0,float *q0,float *v,float *v0,int m,\ int n,float *e,float *f,int k,float *g,float *b,float *dd){extern FILE *file4;FILE *fp;int i,j,l;int pos1,pos2;float a1,a2,d1,d;if(file4==NULL){fp=stdout; /* 输出到屏幕*/}else{fp=file4; /* 输出到文件*/}l=n-1;if(k==1){fprintf(fp,"\n CHANGE OF P0,V**2,P0(I),Q0(I),V0(I) ");fprintf(fp,"\n I P0(I) Q0(I)");}for(i=1;i<=l;i++){a1=0;a2=0;pos1=f1(i);for(j=1;j<=n;j++){/* a1,a2对应课本p185式(4-86)中括号内的式子*/pos2=f2(i,j,n);a1+=g[pos2]*e[f1(j)]-b[pos2]*f[f1(j)];a2+=g[pos2]*f[f1(j)]+b[pos2]*e[f1(j)];}/* 计算式(4-86)(4-87)中的deltaPi */p0[pos1]=p[pos1]-e[pos1]*a1-f[pos1]*a2;if(i <= m){ /* 计算PQ结点中的deltaQi */q0[pos1]=q[pos1]-f[pos1]*a1+e[pos1]*a2;}{ /* 计算PV结点中的deltaVi平方*/v0[pos1]=v[pos1]*v[pos1]-e[pos1]*e[pos1]-f[pos1]*f[pos1];}/* 输出结果*/if(k==1){if(i<m){fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]);}else if(i==m){fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]);fprintf(fp,"\n I P0(I) V0(I)");}else{fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],v0[pos1]);}}}/* 找到deltaP和deltaQ中的最小者，作为收敛指标, 存在dd中*/d=0;for(i=1;i<=l;i++){pos1=f1(i);d1=p0[pos1]>0 ? p0[pos1] : -p0[pos1];if(d<d1)d=d1;}if(i<=m){d1=q0[pos1]>0?q0[pos1]:-q0[pos1];if(d<d1){d=d1;}}}(*dd)=d;}/**************************************************** 本子程序根据节点导纳及电压求Jacoby矩阵,用于求* * 电压修正量,如打印参数K=1,则输出Jacoby矩阵. * * 对应于课本P186式(4-89)(4-90) ****************************************************/void jmcc(int m,int n,int n0,float *e,float *f,float *g,float *b,float *jm,int k){extern FILE *file4;FILE *fp;int i,j,i1,io,i0,ns;int pos1,pos2;if(file4==NULL){fp=stdout;}{fp=file4;}/* 初始化矩阵jm */for(i=1;i<=n0;i++){for(j=1;j<=n0;j++){jm[f2(i,j,n0)]=0;}}ns=n-1; /* 去掉一个平衡结点*//* 计算式(4-89)(4-90) */for(i=1;i<=ns;i++){/* 计算式(4-90) */for(i1=1;i1<=n;i1++){/* pos1是式(4-90)中的j */pos1=f1(i1);/* pos2是式(4-90)中的ij */pos2=f2(i,i1,n);if(i<=m) /* i是PQ结点*/{/* 计算式(4-90)中的Jii等式右侧第一部分*/jm[f2(2*i-1,2*i-1,n0)]+=g[pos2]*f[pos1]+b[pos2]*e[pos1];/* 计算式(4-90)中的Lii等式右侧第一部分*/jm[f2(2*i-1,2*i,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1]; }/* 计算式(4-90)中的Hii等式右侧第一部分*/jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1]; /* 计算式(4-90)中的Nii等式右侧第一部分*/jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]-b[pos2]*e[pos1];}/* pos2是式(4-90)中的ii */pos2=f2(i,i,n);/* pos1是式(4-90)中的i */pos1=f1(i);if(i<=m) /* i是PQ结点*/{/* 计算式(4-90)中的Jii */jm[f2(2*i-1,2*i-1,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1]; /* 计算式(4-90)中的Lii */jm[f2(2*i-1,2*i,n0)]+=g[pos2]*e[pos1]+b[pos2]*f[pos1]; }/* 计算式(4-90)中的Hii */jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]-b[pos2]*f[pos1];/* 计算式(4-90)中的Jii */jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1];if(i>m) /* PV结点*/{/* 计算式(4-90)中的Rii */jm[f2(2*i-1,2*i-1,n0)]=-2*e[pos1];/* 计算式(4-90)中的Sii */jm[f2(2*i-1,2*i,n0)]=-2*f[pos1];}/* 计算式(4-89) */for(j=1;j<=ns;j++){if(j!=i){/* pos1是式(4-89)中的i */pos1=f1(i);/* pos2是式(4-89)中的ij */pos2=f2(i,j,n);/* 计算式(4-89)中的Nij */jm[f2(2*i,2*j,n0)]=b[pos2]*e[pos1]-g[pos2]*f[pos1]; /* 计算式(4-89)中的Hij */jm[f2(2*i,2*j-1,n0)]=-g[pos2]*e[pos1]-b[pos2]*f[pos1]; if(i<=m) /* i是PQ结点*/{/* 计算式(4-89)中的Lij (=-Hij) */jm[f2(2*i-1,2*j,n0)]=-jm[f2(2*i,2*j-1,n0)];/* 计算式(4-89)中的Jij (=Nij) */jm[f2(2*i-1,2*j-1,n0)]=jm[f2(2*i,2*j,n0)];}else /* i是PV结点*/{/* 计算式(4-89)中的Rij (=0) */jm[f2(2*i-1,2*j-1,n0)]=0;/* 计算式(4-89)中的Sij (=0) */jm[f2(2*i-1,2*j,n0)]=0;}}}}if(k!=1){return;}/* 输出Jacoby矩阵*/fprintf(fp,"\n J MATRIX(C)"); for(io=1;io<=n0;io+=5){i1=(io+4)>n0?n0:(io+4);fprintf(fp,"\n");for(j=io;j<=i1;j++){fprintf(fp,"%10d",j);}for(i=1;i<=n0;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i1;j++){fprintf(fp,"%12.6f",jm[f2(i,j,n0)]);}}}fprintf(fp,"\n");}/*********************************************** 本子程序用选列主元素的高斯消元法求解组* * 性方程组求各结点电压修正量,如打印参数K=1,则* * 输出增广矩阵变换中的上三角及电压修正量.如果* * 无唯一解,则给出信息,并停止程序运行. ***********************************************/void sevc ( float a[], int n0, int k, int n1){extern FILE *file4;FILE *fp;int i,j,l,n2,n3,n4,i0,io,j1,i1;float t0,t,c;if(file4==NULL) fp=stdout;else fp=file4;for(i=1;i<=n0;i++){l=i;for(j=i;j<=n0;j++){if( fabs(a[f2(j,i,n1)]) > fabs(a[f2(l,i,n1)]) ){l=j; /* 找到这行中的最大元*/}}if(l!=i){ /* 行列交换*/for (j=i;j<=n1;j++){t=a[f2(i,j,n1)];a[f2(i,j,n1)]=a[f2(l,j,n1)];a[f2(l,j,n1)]=t;}}if (fabs(a[f2(i,i,n1)]-0)<1e-10){ /* 对角元近似于0, 无解*/ printf("\nNo Solution\n");exit (1);}t0=a[f2(i,i,n1)];for(j=i;j<=n1;j++){/* 除对角元*/a[f2(i,j,n1)]/=t0;}if(i==n0){ /* 最后一行，不用消元*/ continue;}/* 消元*/j1=i+1;for(i1=j1;i1<=n0;i1++){c=a[f2(i1,i,n1)];for(j=i;j<=n1;j++){a[f2(i1,j,n1)] -= a[f2(i,j,n1)] *c;}}}if(k==1){ /* 输出上三角矩阵*/fprintf(fp,"\nTrianglar Angmentex Matrix ");for(io=1;io<=n1;io+=5){i0=(io+4)>n1?n1:(io+4);fprintf(fp,"\n");fprintf(fp," ");for(i=io;i<=i0;i++){fprintf(fp,"%12d",i);}for(i=1;i<=n0;i++){fprintf(fp,"\n");fprintf(fp,"%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%15.6f", a[f2(i,j,n1)]);}}}}/* 回代求方程解*/n2=n1-2;for(i=1;i<=n2;i++){n3=n1-i;for(i1=n3;i1<=n0;i1++){n4=n0-i;a[f2(n4,n1,n1)] -= a[f2(i1,n1,n1)]*a[f2(n4,i1,n1)];}}if(k!=1){return;}/* 输出电压修正值*/fprintf(fp,"\nVoltage correction E(i), F(i) :");for(io=1;io<=n0;io+=4){i1=(io+1)/2;i0=((io+3)/2)>(n0/2)?(n0/2):((io+3)/2);fprintf(fp,"\n");for(j=i1;j<=i0;j++){fprintf(fp,"%16d%16d",j,j);}i1 = 2*i0;fprintf(fp,"\n");for(i=io;i<=i1;i++){fprintf(fp,"%15.6f", a[f2(i,n1,n1)]);}}}#define Pi 3.1415927/180void plsc(int n,int l,int m,float g[],float b[],float e[],float f[],\int e1[],int s1[],float g1[],float b1[],float c1[],float c[],\float co[],float p1[],float q1[],float p2[],float q2[],float p3[],\float q3[],float p[],float q[],float v[],float angle[],int k1) {extern FILE *file4;FILE *fp;float t1,t2,st,en,cm,x,y,z,x1,x2,y1,y2;int i,i1,j,m1,ns,pos1,pos2,km;ns=n-1;if(file4==NULL){fp=stdout;}else{fp=file4;}fprintf(fp,"\nTHE RESULT ARE:");if(k1==1){for(i=0;i<n;i++){angle[i]*=Pi;e[i]=v[i]*cos(angle[i]);f[i]=v[i]*sin(angle[i]);}}t1=0.0;t2=0.0;for(i=1;i<=n;i++){pos1=f1(i);pos2=f2(n,i,n);t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1]; }pos1=f1(n);p[pos1]=t1*e[pos1];q[pos1]=-t2*e[pos1];m1=m+1;for(i1=m1;i1<=ns;i1++){t1=0;t2=0;for(i=1;i<=n;i++){pos1=f1(i);pos2=f2(i1,i,n);t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];}pos1=f1(i1);q[pos1]=f[pos1]*t1-e[pos1]*t2;}for(i=0;i<n; i++){cm=co[i];if(cm!=0){q[i]-=(e[i]*e[i]+f[i]*f[i])*cm;}}fprintf(fp,"\nBUS DATA");fprintf(fp,"\nBUS VOLTAGE ANGLE(DEGS.) B US P BUS Q");for(i=0;i<n;i++){v[i]=sqrt(e[i]*e[i]+f[i]*f[i]);x=e[i];y=f[i];z=y/x;angle[i]=atan(z);angle[i]/=Pi;fprintf(fp,"\n%3d%13.5e%15.5f%15.5e%15.5e",i+1,v[i],angle[i],p[ i],q[i]);}fprintf(fp,"\n LINE FLOW ");for(i=1;i<=l;i++){pos1=f1(i);st=s1[pos1];en=e1[pos1];x1=e[f1(st)]*e[f1(st)]+f[f1(st)]*f[f1(st)];x2=e[f1(en)]*e[f1(en)]+f[f1(en)]*f[f1(en)];y1=e[f1(st)]*e[f1(en)]+f[f1(st)]*f[f1(en)];y2=f[f1(st)]*e[f1(en)]-e[f1(st)]*f[f1(en)];p1[pos1]=(x1-y1)*g1[pos1]-y2*b1[pos1];q1[pos1]=-x1*(c1[pos1]+b1[pos1])+y1*b1[pos1]-y2*g1[pos1];p2[pos1]=(x2-y1)*g1[pos1]+y2*b1[pos1];q2[pos1]=-x2*(c1[pos1]+b1[pos1])+y1*b1[pos1]+y2*g1[pos1];for(j=1;j<=n;j++){cm=c[f2(j,i,l)];if(cm!=0.0){km=1;if(en==j){km=2;}if(km==1){q1[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;}else{q2[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;}}}p3[pos1]=p1[pos1]+p2[pos1] ;q3[pos1]=q1[pos1]+q2[pos1] ;fprintf(fp,"\n%2d%8d%11d%13.6e%13.6e%13.6e%13.6e%17d% 11d%13.6e%13.6e",\i,s1[pos1],e1[pos1],p1[pos1],q1[pos1],p3[pos1],q3[pos1],\e1[pos1],s1[pos1],p2[pos1],q2[pos1]);}}【最新整理，下载后即可编辑】。

潮流计算程序(3.0版)说明1特点∙既可以用于高压输电网的潮流解算，又可以用于低压配电网的潮流解算，还可以同时解算输电网加配电网的混合潮流问题。

∙可以处理多平衡（机）节点问题，使用者只需输入各平衡节点的电压幅值和相位角，计算出的系统不平衡功率部分将自动在各平衡节点间进行分摊。

∙可以同时解算相互解列的几个子系统的潮流问题（只要每个子系统内均含有平衡节点），甚至某些子系统退化成孤立节点也不会影响其它子系统的潮流解算。

∙可以同时解算多配电馈线的潮流问题，而无需一个馈线一个馈线地分别计算。

∙由于该程序能够处理多平衡节点问题，故在解算具有多根节点和环状配电网潮流时，无需象以往方法那样只保留一个平衡节点，而将其余平衡节点全部人为地改成PV节点。

∙该版程序为PQ分解法。

2一般变量说明∙LINEMAX 程序所能处理的最大线路数，可以在#define说明语句中进行修改∙GENERA TORMAX 程序所能处理的最大发电机节点数，可以在#define说明语句中进行修改∙LOADMAX 程序所能处理的最大负荷节点数，可以在#define说明语句中进行修改∙NODEMAX 程序所能处理的最大系统节点总数，可以在#define说明语句中进行修改∙SWINGMAX 程序所能处理的最大平衡节点数，可以在#define说明语句中进行修改∙PVMAX 程序所能处理的最大PV节点数，可以在#define说明语句中进行修改∙NODEFACTOR 导纳矩阵的上三角阵中程序所能处理的最大非零非对角元素的个数相对于最大节点数（NODEMAX）的倍数∙Deg_to_Rad 度到弧度的转换系数，在#define说明语句中定义∙Rad_to_Deg 弧度到度的转换系数，在#define说明语句中定义∙Num_Line 系统的实际总线路数∙Num_Gen 系统的实际总发电机（节点）数∙Num_Load 系统的实际总负荷（节点）数∙Num_Node 系统的实际总节点数∙Num_Swing 系统的实际总平衡节点数∙Num_GPV 系统发电机节点中的PV节点数∙Num_GPQ 系统发电机节点中的PQ节点数∙Eps 节点功率失配量值收敛限值∙Iter_Max 迭代次数限值（在不满足电压收敛误差限值的条件下强行中止收敛的最大迭代次数）∙V olIni_Flag 是否读取电压初值标志：0-不读（以0-1启动）；1-读初值∙V olRes_Flag 是否将电压收敛值保存以备以后计算时当启动初值：0-不保存；1-保存3参数结构定义说明∙线路参数结构定义struct Line{int Node_No[2];int Flag;double RXBK[3];}LLine[LINEMAX]其中，Node_No[2]代表线路两端节点号：0—左节点号；1—右节点号。

void factor();/*求因子表*/
void solve(float **B,float *X,int N);/*解方程组*/
void PrtNode();/*打印输出节点参数*/
void ErrorMsg(int Flag);/*错误信息*/
int Node;/*节点数*/
int num;/*原始节点序号*/
kp=0;
for(i=0;i<NP;i++)
{
dPi=dP+i;
Yi=*(Y+i)-i;
Dltai=*(Dlta+i);
*dPi=0;
for(j=0;j<Node;j++)
{
temp=Dltai-*(Dlta+j);
if(i>j)*dPi+=*(V+j)*(Pji);
tP=*(V+j)*(Pij);
tP=*(V+i)*Yij.G-tP;
tP*=*(V+i);
tQ=*(V+j)*(Qij);
tQ-=*(V+i)*(Yij.B-Yij.B0);
tQ*=*(V+i);
}
fprintf(out,"S[%d,%d]=(%10.6f,%10.6f)\n",k+1,m+1,-tP,-tQ)
*(num+i)=k;
fscanf(in,"%d",&k);
}
if(NQ+j!=Node)ErrorMsg(4);
fprintf(out,"【节点参数表】\n");

//////////////////////////////////////////////////////////////////////// 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语言可以通过定义和存储这些参数，并在计算过程中引用它们来实现对这些条件的考虑。

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

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

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

潮流计算，电力学名词，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。

潮流计算是根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。

通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。

待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。

潮流计算程序，相关的原始数据数据数据输入格式如下：%B1是支路参数矩阵，第一列和第二列是起始节点编号和终点节点编号%第三列、第四列、第五列、第六列、第七列、第八列分别为：支路电阻、电抗、电导、电纳、变压器变比、是否有变压器（1为有、0为无）。

%B2为节点参数矩阵，其中第一列到第六列为节点编号；为节点类型；注入有功、注入无功、电压幅值、电压相位。

%“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。

n=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');Y=zeros(n);N=1;%建立节点导纳矩阵for i=1:n1if B1(i,8)==0p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)+B1(i,4))-B1(i,5);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4))+0.5*B1(i,6)+B1(i,5);Y(q,q)=Y(q,q)+1/(B1(i,3)+B1(i,4))+0.5*B1(i,6)+B1(i,5);elsep=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/((B1(i,3)+B1(i,4))*B1(i,7)) -B1(i,5);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4))+B1(i,5);Y(q,q)=Y(q,q)+1/(B1(i,7)^2*(B1(i,3)+B1(i,4)))+B1(i,5);endendYG=real(Y);B=imag(Y);PriS=zeros(2*n-2,1);ImbS=zeros(2*n-2,1);%创建PriS，用于存储初始功率参数h=0;j=0;for i=1:nif i~=isb&B2(i,2)==1h=h+1;for j=1:nPriS(2*h-1,1)=PriS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i, j)*B2(j,5));PriS(2*h,1)=PriS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B 2(j,5));endendendfor i=1:nif i~=isb&B2(i,2)==2h=h+1;for j=1:nPriS(2*h-1,1)=PriS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5));PriS(2*h,1)=PriS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5))endendendPriSU3=zeros(n-h-1,2);%U3存储PV节点的初始电压t=0;for i=1:nif B2(i,2)==2t=t+1;U3(t,1)=B2(i,5);U3(t,2)=B2(i,6);endendU3%ImbS于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:nif i~=isb&B2(i,2)==1h=h+1;ImbS(2*h-1,1)=B2(i,3)-PriS(2*h-1,1);ImbS(2*h,1)=B2(i,4)-PriS(2*h,1);endendt=0;for i=1:nif i~=isb&B2(i,2)==2h=h+1;t=t+1;ImbS(2*h-1,1)=B2(i,3)-PriS(2*h-1,1);ImbS(2*h,1)=U3(t,1)^2+U3(t,2)^2-B2(i,5)^2-B2(i,6)^2;endendImbSI=zeros(n-1,1);%I，存储节点电流参数h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(PriS(2*h-1,1)-PriS(2*h,1)*sqrt(-1))/conj(B2(i,5)+B2(i,6)*sqrt(-1));endendIJacbi=zeros(2*n-2);%Jacbi(雅可比矩阵)h=0;k=0;for i=1:nif B2(i,2)==1h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));elseJacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,2)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=2*B2(i,6);Jacbi(2*h,2*k)=2*B2(i,5);elseJacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendendendJacbi%求解修正方程，获取节点电压的不平衡量ImbU=zeros(2*n-2,1);ImbU=inv(Jacbi)*ImbS;ImbU%修正节点电压j=0;for i=1:nif B2(i,2)==1j=j+1;B2(i,5)=B2(i,5)+ImbU(2*j,1);B2(i,6)=B2(i,6)+ImbU(2*j-1,1);endendfor i=1:nif B2(i,2)==2j=j+1;B2(i,5)=B2(i,5)+ImbU(2*j,1);B2(i,6)=B2(i,6)+ImbU(2*j-1,1);endendB2while abs(max(ImbU))>prPriS=zeros(2*n-2,1);h=0;j=0;for i=1:nif i~=isb&B2(i,2)==1h=h+1;for j=1:nPriS(2*h-1,1)=PriS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i, j)*B2(j,5));PriS(2*h,1)=PriS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B 2(j,5));endendendfor i=1:nif i~=isb&B2(i,2)==2h=h+1;for j=1:nPriS(2*h-1,1)=PriS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5));PriS(2*h,1)=PriS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5))endendendPriS%创建ImbSh=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,2)==1h=h+1;ImbS(2*h-1,1)=B2(i,3)-PriS(2*h-1,1);ImbS(2*h,1)=B2(i,4)-PriS(2*h,1);endendt=0;for i=1:n %对PV节点的处理if i~=isb&B2(i,2)==2h=h+1;t=t+1;ImbS(2*h-1,1)=B2(i,3)-PriS(2*h-1,1);ImbS(2*h,1)=U3(t,1)^2+U3(t,2)^2-B2(i,5)^2-B2(i,6)^2;endendImbS%创建II=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(PriS(2*h-1,1)-PriS(2*h,1)*sqrt(-1))/conj(B2(i,5)+B2(i,6)*sqrt(-1));endendI%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,2)==1h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));elseJacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,2)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=2*B2(i,6);Jacbi(2*h,2*k)=2*B2(i,5);elseJacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendendendJacbiImbU=zeros(2*n-2,1);ImbU=inv(Jacbi)*ImbS;ImbU%修正节点电压j=0;for i=1:nif B2(i,2)==1j=j+1;B2(i,5)=B2(i,5)+ImbU(2*j,1);B2(i,6)=B2(i,6)+ImbU(2*j-1,1);endendfor i=1:nif B2(i,2)==2j=j+1;B2(i,5)=B2(i,5)+ImbU(2*j,1);B2(i,6)=B2(i,6)+ImbU(2*j-1,1);endendB2N=N+1; %迭代次数加1endN。

/*此程序为用牛顿拉夫逊法的直角坐标形式求电力线路的潮流计算程序*//*此程序为一通用程序,按给定的题目输入数据到input文件中，就可求出相关节点的电压相量，支路功率，网络损耗，输电效率平衡节点的功率，pv节点的无功*//*考虑了线路对地导纳支路、平行支路、非标准变比变压器支路、还有无功越限的处理等4种较为复杂的情况*//*本程序通用性很强*/#include <stdio.h>#include<math.h>#include<stdlib.h>void main(){int bus,line,pQ,pv,blzl,transform;/*分别表示节点数，支路数，pQ节点数，pv节点数，对地支路数，变压器数*/int i,j,k,m,n,circle,num,w,r;/*循环变量*/int type[100],node[100],head[100],end[100],branch[100];/*分别表示:节点类型，节点号，支路的首节点，支路的尾节点，支路号*/int t1[100],t2[100];/*中间数组*/float precise;/*计算精度*/float real[100][100],image[100][100],R[100],X[100];/*节点导纳的实部，节点导纳的虚部，支路上的电阻，支路上的电抗*/float p[100],Q[100],p0[100],Q0[100];/*节点的有功功率，节点的无功功率，各节点的给定有功功率初值，各节点的给定无功功率初值*/float e[100],f[100],e0[100],f0[100];/*分别为各节点电压的实部，虚部，给定电压的初始实部和虚部*/float ub[200],ubp[100],ubQ[100],ubu2[100],a[100],b[100];/*分别表示不平衡量，有功不平衡量，无功不平衡量，电压不平衡量,a[]和b[]是中间数组*/floatH[100][100],N[100][100],B[100][100],J[100][100],L[100][100],R1[100][100],S[100][100];/*雅克比矩阵各元素变量*/float a1[100][100],b1[100][100],du[100];/*中间数组和电压修正量*/floatt,s,E[100][100],Ln[100][100],Un[100][100],Bn[100][100],zgjz[100][200],Eu[100][100];/*求逆矩阵的过程中的各数组*/float dumax,ps,Qs,pline[100][100],Qline[100][100],dpline[100][100],dQline[100][100];/*最大电压修正量，平衡节点的有功和无功，支路有功和无功，支路损耗的有功和无功*/ float dpsum,dQsum,efficiency,p1,p2;/*网络总的有功损耗，网络总的无功损耗，输电效率，网络中各节点注入的正的有功的和，网络中各节点注入负的有功的绝对值和*/ float duidi[100],ratio[100],tt1,tt2;/*对地支路的导纳和变压器变比，还有中间变量tt1和tt2*/int x,y,d[100],flag[100];/*flag用来判断平行支路数*/FILE *fp1,*fp2;/*指针变量*/float ee[100],ff[100];/*中间变量*/float Qmax;/*Qmax为无功上限*//*以上为对程序中用到的变量的定义*/fp1=fopen("input.txt","r");/*打开input.txt文件读入数据*/fp2=fopen("output.txt","w");/*打开out.txt文件写入数据*/fprintf(fp2,"\n 华北电力大学\n\n");fprintf(fp2," 电力系统潮流上机班级:电气化0805 姓名:张同学号:200801000528\n");fprintf(fp2,"\n ***************原始数据***************** \n");fprintf(fp2,"============================================================= ======================\n\n");fscanf(fp1,"%d%d%d%d%d%d%f%f",&bus,&line,&pQ,&pv,&blzl,&transform,&precise,&Qmax );fprintf(fp2," 节点数:%d 支路数:%d pQ节点数:%d pv节点数:%d\n 并联对地之路数:%d 变压器数:%d 要求精度:%f \n",bus,line,pQ,pv,blzl,transform,precise);fprintf(fp2,"------------------------------------------------------------------------------------\n");/*下边开始从input.txt文件读入数据*/for(i=1;i<=bus;i++){fscanf(fp1,"%d%d",&node[i],&type[i]);if(type[i]==1)/*type[i]=1,说明第i节点是pQ节点*/fscanf(fp1,"%f%f%f%f",&p0[i],&Q0[i],&e[i],&f[i]);else if(type[i]==2)/*type[i]=1,说明第i节点是pv节点*/fscanf(fp1,"%f%f%f%f%f",&p0[i],&e[i],&f[i],&e0[i],&f0[i]);else if(type[i]==3)/*type[i]=1,说明第i节点是平衡节点*/fscanf(fp1,"%f%f",&e[i],&f[i]);}for(i=1;i<=bus;i++)/*bus表示节点数*/{if(type[i]==1)fprintf(fp2," 节点%d pQ节点%f %f\n",node[i],p0[i],Q0[i]);else if(type[i]==2)fprintf(fp2," 节点%d pv节点%f %f\n",node[i],p0[i],e0[i]);elsefprintf(fp2," 节点%d 平衡节点%f %f\n",node[i],e[i],f[i]);}fprintf(fp2,"------------------------------------------------------------------------------------\n");for(i=1;i<=line;i++){fscanf(fp1,"%d%d%d%f%f%f%d%f",&branch[i],&head[i],&end[i],&R[i],&X[i],&duidi[i],&flag[ i],&ratio[i]);/*读支路信息*/fprintf(fp2," 支路%d %d %d R[%d]=%fX[%d]=%f\n",branch[i],head[i],end[i],branch[i],R[i],branch[i],X[i]); /*打印出给定题目的支路信息*/}/*head里存放支路的首节点，end里存放支路的尾节点*/fprintf(fp2,"============================================================= ======================\n\n");fprintf(fp2,"\n *********************潮流计算开始***************** \n");fclose(fp1);/*关闭文件input.txt文件*//*下边求节点导纳阵开始*//*下面是考虑平行之路时的支路阻抗计算*/for(i=1;i<=line;i++){if(flag[i]==2){R[i]=R[i]/2;X[i]=X[i]/2;duidi[i]=duidi[i]*2;}}/*考虑平行之路的阻抗结束*/for(i=1;i<=bus;i++)for(j=1;j<=bus;j++){real[i][j]=0;/*real是节点导纳的实部*/image[i][j]=0;/*image是节点导纳的虚部*/}for(i=1;i<=line;i++)/*line是总的支路数*/{real[(head[i])][(end[i])]=(-1)*R[i]/(R[i]*R[i]+X[i]*X[i]);image[(head[i])][(end[i])]=X[i]/(R[i]*R[i]+X[i]*X[i]);real[(end[i])][(head[i])]=real[(head[i])][(end[i])];image[(end[i])][(head[i])]=image[(head[i])][(end[i])];for(j=1;j<=bus;j++)if((j==head[i])||(j==end[i])){real[j][j]+=(-1)*real[(head[i])][(end[i])];image[j][j]+=(-1)*image[(head[i])][(end[i])];}}/*下边考虑对地之路*/for(i=1;i<=line;i++){image[head[i]][head[i]]+=duidi[i];image[end[i]][end[i]]+=duidi[i];}/*对地之路考虑结束*//*下边是考虑变压器的非标准变比的程序*/for(i=1;i<=line;i++){tt1=real[(head[i])][(end[i])];tt2=image[(head[i])][(end[i])];real[(head[i])][(end[i])]=real[(head[i])][(end[i])]/fabs(ratio[i]);image[(head[i])][(end[i])]=image[(head[i])][(end[i])]/fabs(ratio[i]);real[(end[i])][(head[i])]=real[(head[i])][(end[i])];image[(end[i])][(head[i])]=image[(head[i])][(end[i])];if(ratio[i]>0){real[head[i]][head[i]]=real[head[i]][head[i]]-tt1*(1-ratio[i])/(ratio[i]*ratio[i])-real[(head[i])][(end[i ])]+tt1;image[head[i]][head[i]]=image[head[i]][head[i]]-tt2*(1-ratio[i])/(ratio[i]*ratio[i])-image[(head[i])] [(end[i])]+tt2;}if(ratio[i]<0){real[end[i]][end[i]]=real[(end[i])][(end[i])]-tt1*(1+ratio[i])/(ratio[i]*ratio[i])-real[(head[i])][(end[i] )]+tt1;image[end[i]][end[i]]=image[(end[i])][(end[i])]-tt2*(1+ratio[i])/(ratio[i]*ratio[i])-image[(head[i])] [(end[i])]+tt2;}}/*考虑变压器的非标准变比的程序结束*/fprintf(fp2,"========================================================= =================================\n\n");fprintf(fp2," 节点导纳阵为:Y=\n");/*打印节点导纳阵*/for(i=1;i<=bus;i++){for(j=1;j<=bus;j++)fprintf(fp2," %f+j%f",real[i][j],image[i][j]);fprintf(fp2,"\n");}/*节点导纳阵已经求出，并已经打印*/fprintf(fp2,"========================================================= ==================================\n\n");for (i=1;i<=bus;i++){ ee[i]=e[i];ff[i]=f[i];}loop : for (i=1;i<=bus;i++){ e[i]=ee[i];f[i]=ff[i];}fprintf(fp2,"=====================================迭代过程==============================================\n\n");fprintf(fp2,"迭代过程中的电压和功率初值为:\n");for(i=1;i<=bus;i++){if(type[i]==1)fprintf(fp2," 节点%d e=%f f=%f p=%f Q=%f\n",node[i],e[i],f[i],p0[i],Q0[i]);else if(type[i]==2)fprintf(fp2," 节点%d e=%f f=%f p=%f\n",node[i],e[i],f[i],p0[i]);}fprintf(fp2,"------------------------------------------------------------------------------------\n");circle=1;/*circle为迭代的次数*//*下边迭代开始*/do{fprintf(fp2,"第%d 次迭代\n",circle);for(i=1;i<=bus;i++){if(type[i]==1){p[i]=0;Q[i]=0;}else if(type[i]==2)p[i]=0;}/*下边是求修正方程的不平衡列向量*/fprintf(fp2,"\n修正方程的不平衡量为:\n");m=1;n=1;for(i=1;i<=bus;i++){if(type[i]==1){for(j=1;j<=bus;j++){p[i]+=(e[i]*(real[i][j]*e[j]-image[i][j]*f[j])+f[i]*(real[i][j]*f[j]+image[i][j]*e[j]));Q[i]+=(f[i]*(real[i][j]*e[j]-image[i][j]*f[j])-e[i]*(real[i][j]*f[j]+image[i][j]*e[j]));}a[m]=p0[i]-p[i];a[(++m)++]=Q0[i]-Q[i];/*a数组存放pQ节点的不平衡量*/ }if(type[i]==2){for(j=1;j<=bus;j++){p[i]+=(e[i]*(real[i][j]*e[j]-image[i][j]*f[j])+f[i]*(real[i][j]*f[j]+image[i][j]*e[j]));}b[n]=p0[i]-p[i];b[(++n)++]=e0[i]*e0[i]+f0[i]*f0[i]-e[i]*e[i]-f[i]*f[i];}/*b数组存放pv节点的不平衡量*/}for(i=1;i<=2*pQ;i++)ub[i]=a[i];/*ub存放的是修正方程的不平衡列向量*/for(i=1;i<=2*pv;i++)ub[i+2*pQ]=b[i];for(i=1;i<=2*bus-2;i++)fprintf(fp2," ub[%d]=%f",i,ub[i]);/*不平衡量已经求出并已经打印*//*下边开始求雅克比矩阵*/fprintf(fp2,"\n------------------------------------------------------------------------------------");fprintf(fp2,"\n雅克比矩阵为:");for(i=1;i<=bus;i++)for(j=1;j<=bus;j++){a1[i][j]=0;b1[i][j]=0;}for(i=1;i<=bus;i++){for(j=1;j<=bus;j++)if(i!=j){a1[i][i]+=(real[i][j]*e[j]-image[i][j]*f[j]);b1[i][i]+=(real[i][j]*f[j]+image[i][j]*e[j]);}a1[i][i]+=(real[i][i]*e[i]-image[i][i]*f[i]);b1[i][i]+=(real[i][i]*f[i]+image[i][i]*e[i]);}for(i=1;i<=bus;i++)for(j=1;j<=bus;j++){if(type[i]==1){H[i][j]=(-1)*image[i][j]*e[i]+real[i][j]*f[i];N[i][j]=real[i][j]*e[i]+image[i][j]*f[i];J[i][j]=(-1)*N[i][j];L[i][j]=H[i][j];}else if(type[i]==2){H[i][j]=(-1)*image[i][j]*e[i]+real[i][j]*f[i];N[i][j]=real[i][j]*e[i]+image[i][j]*f[i];if(i==j){R1[i][j]=2*f[i];S[i][j]=2*e[i];}else if(i!=j){R1[i][j]=0;S[i][j]=0;}}if(i==j){H[i][j]=H[i][j]+b1[i][j];N[i][j]=N[i][j]+a1[i][j];J[i][j]=J[i][j]+a1[i][j];L[i][j]=L[i][j]-b1[i][j];}}num=1;for(i=1;i<=bus;i++){if(type[i]==1){t1[num]=i;num++;}}num=1;{if(type[j]==2){t2[num]=j;num++;}}for(i=1;i<=pQ;i++)d[i]=t1[i];for(i=pQ+1;i<=bus-1;i++)d[i]=t2[i-pQ];for(i=1;i<=2*pQ;i=i+2)for(j=1;j<=2*bus-2;j=j+2){x=d[(i+1)/2];y=d[(j+1)/2];B[i][j]=H[x][y];B[i][j+1]=N[x][y];B[i+1][j]=J[x][y];B[i+1][j+1]=L[x][y];}for(i=2*pQ+1;i<2*bus-1;i=i+2)for(j=1;j<2*bus-1;j=j+2){x=d[(i+1)/2];y=d[(j+1)/2];B[i][j]=H[x][y];B[i][j+1]=N[x][y];B[i+1][j]=R1[x][y];B[i+1][j+1]=S[x][y];}fprintf(fp2,"\n\n");for(i=1;i<2*bus-1;i++){for(j=1;j<2*bus-1;j++)fprintf(fp2," %f",B[i][j]);fprintf(fp2,"\n");}fprintf(fp2,"\n----------------------------------------------------------------------------");/*下面是用LU分解法求雅克比矩阵的逆矩阵的程序*/w=2*bus-2;/*雅克比矩阵的阶数*/for(j=1;j<=w;j++){if(i==j){E[i][j]=1;}elseE[i][j]=0;}/*单位矩阵已经生成*/for(i=1;i<=w;i++)for(j=1;j<=w;j++)zgjz[i][j]=B[i][j];for(i=1;i<=w;i++)for(j=w+1;j<=2*w;j++)zgjz[i][j]=E[i][j-w];for(i=1;i<=2*w;i++)for(j=1;j<=w;j++)Un[1][i]=zgjz[1][i];for(k=2;k<=w;k++)Ln[k][1]=zgjz[k][1]/Un[1][1];for(r=2;r<=w;r++){for(i=r;i<=2*w;i++){t=0;for(k=1;k<=r-1;k++)t+=Ln[r][k]*Un[k][i];Un[r][i]=zgjz[r][i]-t;}for(j=r+1;j<=w;j++){s=0;for(k=1;k<=r-1;k++)s+=Ln[j][k]*Un[k][r];Ln[j][r]=((zgjz[j][r]-s)/Un[r][r]);}}for(i=1;i<=w;i++)for(j=1;j<i;j++)Un[i][j]=0;for(i=1;i<=w;i++)for(j=1;j<=w;j++)Eu[i][j]=Un[i][j+w];for(j=w;j>=1;j--)Bn[w][j]=Eu[w][j]/Un[w][w];for(i=w-1;i>=1;i--)for(j=1;j<=w;j++){t=0;for(k=i+1;k<=w;k++)t+=Un[i][k]*Bn[k][j];for(k=i+1;k<=w;k++)Bn[i][j]=(Eu[i][j]-t)/Un[i][i];/*Bn[i][j]就是原雅克比矩阵的逆矩阵*/ }/*求逆矩阵的程序结束*/for(i=1;i<=w;i++)du[i]=0;for(i=1;i<=w;i++)for(j=1;j<=w;j++)du[i]+=Bn[i][j]*ub[j];fprintf(fp2,"\n修正量为:\n");for(i=1;i<=w;i++)fprintf(fp2," du[%d]=%.6f",i,du[i]);fprintf(fp2,"\n-------------------------------------------------------------------------\n");for(i=1;i<=pQ;i++){f[d[i]]+=du[(2*i-1)];e[d[i]]+=du[(2*i)];fprintf(fp2," f[%d]=%f e[%d]=%f\n",d[i],f[d[i]],d[i],e[d[i]]);}for(i=pQ+1;i<=bus-1;i++){f[d[i]]+=du[2*i-1];e[d[i]]+=du[2*i];fprintf(fp2," f[%d]=%f e[%d]=%f\n",d[i],f[d[i]],d[i],e[d[i]]);}dumax=fabs(du[1]);for(i=2;i<=w;i++){if(fabs(du[i])>dumax)dumax=fabs(du[i]);}circle++;}while(dumax>=0.00010 && circle<=50);fprintf(fp2,"\n-------------------------------------------------------------------------");fprintf(fp2,"\n ***************迭代过程结束*****************\n");fprintf(fp2,"============================================================= ===========================");fprintf(fp2,"\n平衡节点的功率为:\n");ps=0;Qs=0;for(i=1;i<=bus;i++){if(type[i]==3){for(j=1;j<=bus;j++){ps+=(real[i][j]*e[j]-image[i][j]*f[j]);Qs+=((-1)*real[i][j]*f[j]-image[i][j]*e[j]);}ps=ps*e[i];Qs=Qs*e[i];p0[i]=ps;Q0[i]=Qs;fprintf(fp2," S=%f+j%f\n",ps,Qs);}}fprintf(fp2,"pv节点的无功功率为:\n");for(i=1;i<=bus;i++){if(type[i]==2){ Q0[i]=0;for(j=1;j<=bus;j++)Q0[i]+=(f[i]*(real[i][j]*e[j]-image[i][j]*f[j])-e[i]*(real[i][j]*f[j]+image[i][j]*e[j]));fprintf(fp2," Q0[%d]=%f\n",i,Q0[i]);}}/*下边是考虑无功越限情况下的程序*/for(i=1;i<=bus;i++){if(type[i]==2 && Q0[i]>Qmax){ fprintf(fp2,"\n************************无功越限情况下的重新迭代过程*********************\n");pQ++;pv--;Q0[i]=Qmax;/*如果越限，则给pv节点的无功功率赋值为给定的Qmax*/type[i]=1;for(j=1;j<=bus;j++){p[i]=0;Q[i]=0;}goto loop;/*如果pv节点的无功功率超出了给定的Qmax，那么就返回到loop 重新进行迭代过程*/}}/*无功越限程序考虑完毕*//*下边是求线路上的功率的一段程序*/fprintf(fp2,"线路上的功率为:\n");for(i=1;i<=line;i++){real[head[i]][0]=0;real[end[i]][0]=0;image[end[i]][0]=duidi[i];image[head[i]][0]=duidi[i];}for(i=1;i<=line;i++){if(flag[i]==1)/*这是考虑单回线运行时支路的功率和损耗情况*/{pline[head[i]][end[i]]=((real[head[i]][0]+real[head[i]][end[i]])*(e[head[i]]*e[head[i]]+f[head[i]]*f [head[i]])-(e[head[i]]*e[end[i]]+f[head[i]]*f[end[i]])*real[head[i]][end[i]]-(f[head[i]]*e[end[i]]-e[ head[i]]*f[end[i]])*image[head[i]][end[i]])*(-1);Qline[head[i]][end[i]]=((-1)*(image[head[i]][0]+image[head[i]][end[i]])*(e[head[i]]*e[head[i]]+f[ head[i]]*f[head[i]])+(e[head[i]]*e[end[i]]+f[head[i]]*f[end[i]])*image[head[i]][end[i]]+(e[head[i] ]*f[end[i]]-f[head[i]]*e[end[i]])*real[head[i]][end[i]])*(-1);pline[end[i]][head[i]]=((real[end[i]][0]+real[end[i]][head[i]])*(e[end[i]]*e[end[i]]+f[end[i]]*f[end [i]])-(e[head[i]]*e[end[i]]+f[head[i]]*f[end[i]])*real[end[i]][head[i]]-(f[end[i]]*e[head[i]]-e[end[i] ]*f[head[i]])*image[end[i]][head[i]])*(-1);Qline[end[i]][head[i]]=((-1)*(image[end[i]][0]+image[end[i]][head[i]])*(e[end[i]]*e[end[i]]+f[en d[i]]*f[end[i]])+(e[head[i]]*e[end[i]]+f[head[i]]*f[end[i]])*image[end[i]][head[i]]-(f[end[i]]*e[he ad[i]]+e[end[i]]*f[head[i]])*real[end[i]][head[i]])*(-1);dpline[head[i]][end[i]]=pline[head[i]][end[i]]+pline[end[i]][head[i]];dQline[head[i]][end[i]]=Qline[head[i]][end[i]]+Qline[head[i]][end[i]];fprintf(fp2," S[%d][%d]=%f +j %f\n",head[i],end[i],pline[head[i]][end[i]],Qline[head[i]][end[i]]);fprintf(fp2," S[%d][%d]=%f +j %f\n",end[i],head[i],pline[end[i]][head[i]],Qline[end[i]][head[i]]);}else if(flag[i]==2)/*这是考虑双回线运行时支路的功率和损耗情况*/{pline[head[i]][end[i]]=pline[head[i]][end[i]]/2;Qline[head[i]][end[i]]=Qline[head[i]][end[i]]/2;pline[end[i]][head[i]]=pline[end[i]][head[i]]/2;Qline[end[i]][head[i]]=Qline[end[i]][head[i]]/2;dpline[head[i]][end[i]]=dpline[head[i]][end[i]]/2;dQline[head[i]][end[i]]=dQline[head[i]][end[i]]/2;fprintf(fp2," S[%d][%d]=%f +j %f\n",head[i],end[i],pline[head[i]][end[i]],Qline[head[i]][end[i]]);fprintf(fp2," S[%d][%d]=%f +j %f\n",end[i],head[i],pline[end[i]][head[i]],Qline[end[i]][head[i]]);}}fprintf(fp2,"\n-------------------------------------------------------------------------------------------\n");fprintf(fp2,"考虑平行之路和接地之路情况时\n单条线路损耗的功率为:\n");for(i=1;i<=line;i++)fprintf(fp2,"\n 线路%d 上的功率损耗为: %f +j %f",i,dpline[head[i]][end[i]],dQline[head[i]][end[i]]);fprintf(fp2,"\n-------------------------------------------------------------------------------------------\n");/*下边是求网络总损耗的程序*/fprintf(fp2,"网络总损耗为:\n");dpsum=0;dQsum=0;for(i=1;i<=bus;i++)dpsum+=p0[i];for(i=1;i<=bus;i++)dQsum+=Q0[i];fprintf(fp2," ds= %f +j %f\n",dpsum,dQsum);/*网络总损耗完毕*/fprintf(fp2,"-----------------------------------------------------------------------------------------------\n");/*下边是求数点效率的程序*/fprintf(fp2,"输电效率为:\n");p1=0;p2=0;for(i=1;i<=bus;i++){if(type[i]!=3){if (p0[i]>0)p1+=p0[i];elsep2+=(-1)*p0[i];}}p1=p1+ps;efficiency=p2/p1*100;/*输电效率完毕*/fprintf(fp2,"efficiency=%.4f％",efficiency);/*efficiency是输电效率*/fprintf(fp2,"\n============================================================ ===================================\n");fprintf(fp2," **********************潮流计算结束*********************");fclose(fp2);/*关闭output.txt文件*/printf("!!******计算结果存放于out.txt文件中******!!\n");}。

电力系统潮流计算注：这是一个基于N－R法的潮流计算通用程序，仅提供了子程序，需要做些处理才能成为一个可运行的计算程序！此程序非我原创，仅与大家共享！！！/******************************************************************* * 这里提供的是电力系统潮流计算机解法的五个子程序,采用的方法是 ** Newton_Raphson法.** 程序中所用的变量说明如下:** N:网络节点总数. M:网络的PQ节点数. ** L:网络的支路总数. N0:雅可比矩阵的行数. ** N1:N0+1K:打印开关.K=1,则打印;否则,不打印.** K1:子程序PLSC中判断输入电压的形式.K1=1,则为极座标形式.否则** 为直角坐标形式.** D:有功及无功功率误差的最大值. * * G(I,J):Ybus的电导元素(实部).** B(I,J):Ybus的电纳元素(虚部).** G1(I) :第I支路的串联电导. B1(I):第I支路的串联电纳. ** C1(I) :第I支路的pie型对称接地电纳. ** C(I,J):第I节点J支路不对称接地电纳. ** CO(I) :第I节点的接地电纳.** S1(I) :第I节点的起始节点号. E1(I):第I节点的终止节点号. ** P(I) :第I节点的注入有功功率. Q(I):第I节点的注入无功功率.** P0(I) :第I节点有功功率误差. Q0(I):第I节点无功功率误差. ** V0(I) :第I节点(PV节点)的电压误差(平方误差). ** V(I) :第I节点的电压误差幅值. * * E(I) :第I节点的电压的实部. F(I):第I节点的电压的虚部. ** JM(I,J):Jacoby矩阵的第I行J列元素. ** A(I,J):修正方程的增广矩阵,三角化矩阵的第I行J列元素,运算结 ** 束后A矩阵的最后一列存放修正的解. ** P1(I) :第I支路由S1(I)节点注入的有功功率. ** Q1(I) :第I支路由S1(I)节点注入的无功功率. ** P2(I) :第I支路由E1(I)节点注入的有功功率. ** Q2(I) :第I支路由E1(I)节点注入的无功功率. ** P3(I) :第I支路的有功功率损耗. * * Q3(I) :第I支路的无功功率损耗. * * ANGLE(I):第I节点电压的角度.********************************************************************/ #include <math.h>#include <stdio.h>#define f1(i) (i-1)/* 把习惯的一阶矩阵的下标转化为C语言数组下标*/#define f2(i,j,n) ((i-1)*(n)+j-1)/* 把习惯的二阶矩阵的下标转化为C语言数组下标*//***************************************************** 本子程序根据所给的支路导纳及有关信息,形成结点 ** 导纳矩阵,如打印参数K=1,则输出电导矩阵G和电纳矩B *****************************************************/void ybus(int n,int l,int m,float *g,float *b,float *g1,float *b1,float *c1,\float *c,float *co,int k,int *s1,int *e1) {extern FILE *file4;FILE *fp;int i,j,io,i0;int pos1,pos2;int st,en;if(file4==NULL){fp=stdout;}else{fp=file4; /* 输出到文件 */}/* 初始化矩阵G,B */for(i=1;i<=n;i++){for(j=1;j<=n;j++){pos2=f2(i,j,n);g[pos2]=0;b[pos2]=0;}}/* 计算支路导纳 */for(i=1;i<=l;i++){/* 计算对角元 */pos1=f1(i);st=s1[pos1];en=e1[pos1];pos2=f2(st,st,n);g[pos2]+=g1[pos1];b[pos2]+=b1[pos1]+c1[pos1];pos2=f2(en,en,n);g[pos2]+=g1[pos1];b[pos2]+=b1[pos1]+c1[pos1];/* 计算非对角元 */pos2=f2(st,en,n);g[pos2]-=g1[pos1];b[pos2]-=b1[pos1];g[f2(en,st,n)]=g[f2(st,en,n)];b[f2(en,st,n)]=b[f2(st,en,n)];}/* 计算接地支路导纳 */for(i=1;i<=n;i++){/* 对称部分*/b[f2(i,i,n)]+=co[f1(i)];/* 非对称部分*/for(j=1;j<=l;j++){b[f2(i,i,n)]+=c[f2(i,j,l)];}}if(k!=1){return; /* 如果K不为 1,则返回;否则,打印导纳矩阵 */ }fprintf(fp,"\n BUS ADMITTANCE MATRIXY(BUS):");fprintf(fp,"\n ******************* ARRAY G ********************"); for(io=1;io<=n;io+=5){i0=(io+4)>n?n:(io+4);fprintf(fp,"\n");for(j=io;j<=i0;j++){fprintf(fp,"%13d",j);}for(i=1;i<=n;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%13.6f",g[f2(i,j,n)]);}}fprintf(fp,"\n");}fprintf(fp,"\n ******************* ARRAY B ********************");for(io=1;io<=n;io+=5){i0=(io+4)>n?n:(io+4);fprintf(fp,"\n");for(j=io;j<=i0;j++){fprintf(fp,"%13d",j);}for(i=1;i<=n;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%13.6f",b[f2(i,j,n)]);}}fprintf(fp,"\n");}fprintf(fp,"\n************************************************");}/******************************************** 本子程序根据所给的功率及电压等数据 ** 求出功率及电压误差量,并返回最大有功功率 ** 以用于与给定误差比较.如打印参数K=1,则输 ** 出P0,Q0(对PQ结点),V0(对PV结点). ********************************************/void dpqc(float *p,float *q,float *p0,float *q0,float *v,float *v0,int m,\int n,float *e,float *f,int k,float *g,float *b,float *dd){extern FILE *file4;FILE *fp;int i,j,l;int pos1,pos2;float a1,a2,d1,d;if(file4==NULL){fp=stdout; /* 输出到屏幕 */}else{fp=file4; /* 输出到文件*/}l=n-1;if(k==1){fprintf(fp,"\n CHANGE OFP0,V**2,P0(I),Q0(I),V0(I) ");fprintf(fp,"\n I P0(I)Q0(I)");}for(i=1;i<=l;i++){a1=0;a2=0;pos1=f1(i);for(j=1;j<=n;j++){/* a1,a2对应课本p185式(4-86)中括号内的式子 */pos2=f2(i,j,n);a1+=g[pos2]*e[f1(j)]-b[pos2]*f[f1(j)];a2+=g[pos2]*f[f1(j)]+b[pos2]*e[f1(j)];}/* 计算式(4-86)(4-87)中的deltaPi */p0[pos1]=p[pos1]-e[pos1]*a1-f[pos1]*a2;if(i <= m){ /* 计算PQ结点中的deltaQi */q0[pos1]=q[pos1]-f[pos1]*a1+e[pos1]*a2;}else{ /* 计算PV结点中的deltaVi平方*/v0[pos1]=v[pos1]*v[pos1]-e[pos1]*e[pos1]-f[pos1]*f[pos1];}/* 输出结果 */if(k==1){if(i<m){fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]);}else if(i==m){fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]);fprintf(fp,"\n I P0(I)V0(I)");}else{fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],v0[pos1]);}}}/* 找到deltaP和deltaQ中的最小者，作为收敛指标, 存在dd中 */d=0;for(i=1;i<=l;i++){pos1=f1(i);d1=p0[pos1]>0 ? p0[pos1] : -p0[pos1];if(d<d1){d=d1;}if(i<=m){d1=q0[pos1]>0?q0[pos1]:-q0[pos1];if(d<d1){d=d1;}}}(*dd)=d;}/**************************************************** 本子程序根据节点导纳及电压求Jacoby矩阵,用于求** 电压修正量,如打印参数K=1,则输出Jacoby矩阵. ** 对应于课本P186式(4-89)(4-90) ****************************************************/void jmcc(int m,int n,int n0,float *e,float *f,float *g,float *b,float *jm,int k){extern FILE *file4;FILE *fp;int i,j,i1,io,i0,ns;int pos1,pos2;if(file4==NULL){fp=stdout;}else{fp=file4;}/* 初始化矩阵jm */for(i=1;i<=n0;i++){for(j=1;j<=n0;j++){jm[f2(i,j,n0)]=0;}}ns=n-1; /* 去掉一个平衡结点 *//* 计算式(4-89)(4-90) */for(i=1;i<=ns;i++){/* 计算式(4-90) */for(i1=1;i1<=n;i1++){/* pos1是式(4-90)中的j */pos1=f1(i1);/* pos2是式(4-90)中的ij */pos2=f2(i,i1,n);if(i<=m) /* i是PQ结点 */{/* 计算式(4-90)中的Jii等式右侧第一部分 */jm[f2(2*i-1,2*i-1,n0)]+=g[pos2]*f[pos1]+b[pos2]*e[pos1 ];/* 计算式(4-90)中的Lii等式右侧第一部分 */jm[f2(2*i-1,2*i,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1];}/* 计算式(4-90)中的Hii等式右侧第一部分 */jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1];/* 计算式(4-90)中的Nii等式右侧第一部分 */jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]-b[pos2]*e[pos1];}/* pos2是式(4-90)中的ii */pos2=f2(i,i,n);/* pos1是式(4-90)中的i */pos1=f1(i);if(i<=m) /* i是PQ结点 */{/* 计算式(4-90)中的Jii */jm[f2(2*i-1,2*i-1,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1] ;/* 计算式(4-90)中的Lii */jm[f2(2*i-1,2*i,n0)]+=g[pos2]*e[pos1]+b[pos2]*f[po s1];}/* 计算式(4-90)中的Hii */jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]-b[pos2]*f[pos1];/* 计算式(4-90)中的Jii */jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1];if(i>m) /* PV结点 */{/* 计算式(4-90)中的Rii */jm[f2(2*i-1,2*i-1,n0)]=-2*e[pos1];/* 计算式(4-90)中的Sii */jm[f2(2*i-1,2*i,n0)]=-2*f[pos1];}/* 计算式(4-89) */for(j=1;j<=ns;j++){if(j!=i){/* pos1是式(4-89)中的i */pos1=f1(i);/* pos2是式(4-89)中的ij */pos2=f2(i,j,n);/* 计算式(4-89)中的Nij */jm[f2(2*i,2*j,n0)]=b[pos2]*e[pos1]-g[pos2]*f[pos1];/* 计算式(4-89)中的Hij */jm[f2(2*i,2*j-1,n0)]=-g[pos2]*e[pos1]-b[pos2]*f[pos1];if(i<=m) /* i是PQ结点 */{/* 计算式(4-89)中的Lij (=-Hij) */jm[f2(2*i-1,2*j,n0)]=-jm[f2(2*i,2*j-1,n0)];/* 计算式(4-89)中的Jij (=Nij) */jm[f2(2*i-1,2*j-1,n0)]=jm[f2(2*i,2*j,n0)];}else /* i是PV结点 */{/* 计算式(4-89)中的Rij (=0) */jm[f2(2*i-1,2*j-1,n0)]=0;/* 计算式(4-89)中的Sij (=0) */jm[f2(2*i-1,2*j,n0)]=0;}}}}if(k!=1){return;}/* 输出Jacoby矩阵 */fprintf(fp,"\n J MATRIX(C)"); for(io=1;io<=n0;io+=5){i1=(io+4)>n0?n0:(io+4);fprintf(fp,"\n");for(j=io;j<=i1;j++){fprintf(fp,"%10d",j);}for(i=1;i<=n0;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i1;j++){fprintf(fp,"%12.6f",jm[f2(i,j,n0)]);}}}fprintf(fp,"\n");}/*********************************************** 本子程序用选列主元素的高斯消元法求解组 ** 性方程组求各结点电压修正量,如打印参数K=1,则** 输出增广矩阵变换中的上三角及电压修正量.如果** 无唯一解,则给出信息,并停止程序运行. ***********************************************/void sevc ( float a[], int n0, int k, int n1){extern FILE *file4;FILE *fp;int i,j,l,n2,n3,n4,i0,io,j1,i1;float t0,t,c;if(file4==NULL) fp=stdout;else fp=file4;for(i=1;i<=n0;i++){l=i;for(j=i;j<=n0;j++){if( fabs(a[f2(j,i,n1)]) > fabs(a[f2(l,i,n1)]) ){l=j; /* 找到这行中的最大元 */}}if(l!=i){ /* 行列交换 */for (j=i;j<=n1;j++){t=a[f2(i,j,n1)];a[f2(i,j,n1)]=a[f2(l,j,n1)];a[f2(l,j,n1)]=t;}}if (fabs(a[f2(i,i,n1)]-0)<1e-10){ /* 对角元近似于0, 无解 */printf("\nNo Solution\n");exit (1);}t0=a[f2(i,i,n1)];for(j=i;j<=n1;j++){/* 除对角元 */a[f2(i,j,n1)]/=t0;}if(i==n0){ /* 最后一行，不用消元 */continue;}/* 消元 */j1=i+1;for(i1=j1;i1<=n0;i1++){c=a[f2(i1,i,n1)];for(j=i;j<=n1;j++){a[f2(i1,j,n1)] -= a[f2(i,j,n1)] *c;}}}if(k==1){ /* 输出上三角矩阵 */fprintf(fp,"\nTrianglar Angmentex Matrix ");for(io=1;io<=n1;io+=5){i0=(io+4)>n1?n1:(io+4);fprintf(fp,"\n");fprintf(fp," ");for(i=io;i<=i0;i++){fprintf(fp,"%12d",i);}for(i=1;i<=n0;i++){fprintf(fp,"\n");fprintf(fp,"%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%15.6f", a[f2(i,j,n1)]);}}}}/* 回代求方程解 */n2=n1-2;for(i=1;i<=n2;i++){n3=n1-i;for(i1=n3;i1<=n0;i1++){n4=n0-i;a[f2(n4,n1,n1)] -= a[f2(i1,n1,n1)]*a[f2(n4,i1,n1)];}}if(k!=1){return;}/* 输出电压修正值 */fprintf(fp,"\nVoltage correction E(i), F(i) :");for(io=1;io<=n0;io+=4){i1=(io+1)/2;i0=((io+3)/2)>(n0/2)?(n0/2):((io+3)/2);fprintf(fp,"\n");for(j=i1;j<=i0;j++){fprintf(fp,"%16d%16d",j,j);}i1 = 2*i0;fprintf(fp,"\n");for(i=io;i<=i1;i++){fprintf(fp,"%15.6f", a[f2(i,n1,n1)]);}}}#define Pi 3.1415927/180void plsc(int n,int l,int m,float g[],float b[],float e[],float f[],\int e1[],int s1[],float g1[],floatb1[],float c1[],float c[],\float co[],float p1[],float q1[],floatp2[],float q2[],float p3[],\float q3[],float p[],float q[],floatv[],float angle[],int k1){extern FILE *file4;FILE *fp;float t1,t2,st,en,cm,x,y,z,x1,x2,y1,y2;int i,i1,j,m1,ns,pos1,pos2,km;ns=n-1;if(file4==NULL){fp=stdout;}else{fp=file4;}fprintf(fp,"\nTHE RESULT ARE:");if(k1==1){for(i=0;i<n;i++){angle[i]*=Pi;e[i]=v[i]*cos(angle[i]);f[i]=v[i]*sin(angle[i]);}}t1=0.0;t2=0.0;for(i=1;i<=n;i++){pos1=f1(i);pos2=f2(n,i,n);t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];}pos1=f1(n);p[pos1]=t1*e[pos1];q[pos1]=-t2*e[pos1];m1=m+1;for(i1=m1;i1<=ns;i1++){t1=0;t2=0;for(i=1;i<=n;i++){pos1=f1(i);pos2=f2(i1,i,n);t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];}pos1=f1(i1);q[pos1]=f[pos1]*t1-e[pos1]*t2;}for(i=0;i<n; i++){cm=co[i];if(cm!=0){q[i]-=(e[i]*e[i]+f[i]*f[i])*cm;}}fprintf(fp,"\nBUS DATA");fprintf(fp,"\nBUS VOLTAGE ANGLE(DEG S.) BUS P BUS Q");for(i=0;i<n;i++){v[i]=sqrt(e[i]*e[i]+f[i]*f[i]);x=e[i];y=f[i];z=y/x;angle[i]=atan(z);angle[i]/=Pi;fprintf(fp,"\n%3d%13.5e%15.5f%15.5e%15.5e",i+1,v[i],angle[ i],p[i],q[i]);}fprintf(fp,"\n LINE FLOW ");for(i=1;i<=l;i++){pos1=f1(i);st=s1[pos1];en=e1[pos1];x1=e[f1(st)]*e[f1(st)]+f[f1(st)]*f[f1(st)];x2=e[f1(en)]*e[f1(en)]+f[f1(en)]*f[f1(en)];y1=e[f1(st)]*e[f1(en)]+f[f1(st)]*f[f1(en)];y2=f[f1(st)]*e[f1(en)]-e[f1(st)]*f[f1(en)];p1[pos1]=(x1-y1)*g1[pos1]-y2*b1[pos1];q1[pos1]=-x1*(c1[pos1]+b1[pos1])+y1*b1[pos1]-y2*g1[pos1];p2[pos1]=(x2-y1)*g1[pos1]+y2*b1[pos1];q2[pos1]=-x2*(c1[pos1]+b1[pos1])+y1*b1[pos1]+y2*g1[pos1];for(j=1;j<=n;j++){cm=c[f2(j,i,l)];if(cm!=0.0){km=1;if(en==j){km=2;}if(km==1){q1[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;}else{q2[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;}}}p3[pos1]=p1[pos1]+p2[pos1] ;q3[pos1]=q1[pos1]+q2[pos1] ;fprintf(fp,"\n%2d%8d%11d%13.6e%13.6e%13.6e%13.6e%17d%11d%1 3.6e%13.6e",\i,s1[pos1],e1[pos1],p1[pos1],q1[pos1],p3[pos1],q3[pos1 ],\e1[pos1],s1[pos1],p2[pos1],q2[pos1]);}}。

#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<<set w(14)<<Vl<<setw(14)<<low<<endl;fout.close();}。

电力系统潮流及短路电流计算程序以下是一个简单的电力系统潮流计算程序的框架：1.输入数据准备阶段：-输入潮流计算的电力系统拓扑结构，包括各节点之间的连接关系、导线电阻、电抗等信息。

-输入电力系统的负荷信息，包括负荷节点、负荷大小、负荷类型等。

-输入电力系统的发电机信息，包括发电机节点、发电机类型、发电机容量等。

2.潮流计算阶段：-初始条件设置：给定电力系统中各节点的初始电压、相角等信息。

-节点功率方程求解：根据电力系统的拓扑结构和发电机、负荷信息，建立节点功率方程。

-潮流计算迭代：使用牛顿-拉夫逊法等迭代算法求解节点功率方程，得到各节点的电压、相角等参数。

3.潮流计算结果输出阶段：-输出各节点的电压、相角、有功功率、无功功率等参数。

-输出各支路的电流、功率损耗等参数。

-输出系统的功率平衡情况。

4.短路电流计算阶段：-输入短路电流计算的电力系统拓扑结构。

-输入短路电流计算的负荷信息。

-输入短路电流计算的电源信息。

-使用KVL（电压法）或KCL（电流法）等方法计算各节点短路电流。

5.短路电流计算结果输出阶段：-输出各节点的短路电流大小。

-输出各支路的短路电流大小。

以上只是一个电力系统潮流及短路电流计算程序的大致流程框架，具体实现细节和算法选择还需要根据具体情况进行进一步的设计和开发。

在实际应用中，还需要考虑各种特殊情况和计算优化方法，以提高计算速度和准确性。

总之，电力系统潮流及短路电流计算程序是电力工程师在设计和运行电力系统中不可或缺的工具，它能够帮助工程师快速了解系统的运行状态和电流分布情况，以便进行系统优化和安全评估。

用C 语言进行复杂电力系统的潮流计算杨汉明电力系统的潮流计算是对电力系统分析的最基本步骤也是最重要的步骤，是指在一定的系统结构和运行条件下，确定系统运行状态的计算，也即是对各母线（节点）电压，各元件（支路）传输电线或功率的计算。

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

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

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

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

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

对角元Y ii 称为自导纳，它等于与该节点I 直接相连的所有支路导纳总和；非对角元Y ij （i ≠j ）称为互导纳或转移导纳，它等于连结节点I ，j 支路导纳的负数，且有Y ij =Y ji ，当节点I ，j 之间没有支路直接相连时，Y ij =Y ji =0。

#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();}。