5节点电力系统潮流计算程序
电力系统潮流计算完整c语言程序(含网损计算的最终版)
{
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;
电力系统潮流计算
3.2.1 节点电压方程与节点导纳矩阵和阻抗矩阵
将节点电压法应用于电力系统潮流计算，变量为节点电压与节
点注入电流。通常以大地作为电压幅值的参考（U0 = 0），以
系统中某一指定母线的电压角度作为电压相角的参考，以支路
导纳作为电力网的参数进行计算。节点注入电流规定为流向网
络为正，流出为负。
Pmax P
表征年有功负荷曲线特点的两个指标
0
年最大负荷利用小时数 Tmax
t Tmax 8760
根据年负荷曲线，可求得全年所需电能：
8760
A 0
Pdt MWh
定义年最大负荷（最大值 Pmax）利用小时： Tmax
A Pmax
h
Tmax 越大，负荷曲线越平坦
负荷曲线为一水平线时， Tmax 达到最大值8760 (h)
2
1 ZT1
2
Zl
T2
34
3
ZT2 4
YT3
Yl /2
YT2
已知末端功率和电压， 计算网上潮流分布。
1 ZT1 2 Zl
3 ZT2 4
已知始端功率和电压， 计算网上潮流分布。
Y20
Y30
已知末端功率和始端电 压，计算网上的潮流。
不管哪种情况，先作等值电路
3.1.3 辐射形网络的分析计算
1）已知末端功率、电压 利用前面的方法，从末端逐级 往上推算，直至求得各要求的量。
Pm(t)
损耗称年电能损耗，是电网运行经
济性的指标。
Pmi
1）年电能损耗的准确计算方法
已知各负荷的年有功和无功负荷曲线 时，理论上可准确计算年电能损耗。
8760小时分为 n 段，第 i 时段时间为 Dti (h)，全网功率损耗为DPi (MW)，则 全网年电能损耗为
牛顿拉夫逊法计算潮流步骤
牛顿拉夫逊法计算潮流步骤牛顿拉夫逊法（Newton-Raphson method）是一种用于求解非线性方程组的迭代方法，它可以用来计算电力系统潮流的解。
潮流计算是电力系统规划和运行中的重要任务，它的目标是求解电力系统中各节点的电压幅值和相角，以及线路的功率流向等参数，用于分析电力系统的稳定性和安全性，以及进行电力系统规划和调度。
下面是使用牛顿拉夫逊法计算潮流的一般步骤：步骤1：初始化首先，需要对电力系统的各个节点（包括发电机节点和负荷节点）的电压幅值和相角进行初始化，一般可以使用其中一种估计值或者历史数据作为初始值。
步骤2：建立潮流方程根据电力系统的潮流计算模型，可以建立节点电压幅值和相角的平衡方程，一般采用节点注入功率和节点电压的关系来表示。
潮流方程一般是一个非线性方程组，包含了各个节点之间的复杂关系。
步骤3：线性化方程组将潮流方程组进行线性化处理，一般采用泰勒展开的方法，将非线性方程组变为线性方程组。
线性化的过程需要计算雅可比矩阵，即方程组中的系数矩阵。
步骤4：求解线性方程组利用线性方程组的求解方法，比如高斯消元法或LU分解法等，求解线性方程组，得到电压幅值和相角的修正量。
步骤5：更新节点电压根据线性方程组的解，更新各个节点的电压幅值和相角，得到新的节点电压。
步骤6：检查收敛性判断节点电压的修正量是否小于设定的收敛阈值，如果满足收敛条件，则停止迭代；否则，返回步骤3，循环进行线性化方程组和线性方程组的求解。
步骤7：输出结果当潮流计算收敛时，输出最终的节点电压幅值和相角，以及线路的功率流向等参数。
牛顿拉夫逊法是一种高效、快速且收敛性良好的方法，广泛应用于电力系统潮流计算。
在实际应用中，可能会遇到多次迭代或者收敛性不好的情况，此时可以采用退火技术或其他优化算法进行改进。
此外，牛顿拉夫逊法的计算也可以并行化，利用多核处理器或者分布式计算集群来加速计算过程。
总之，牛顿拉夫逊法是一种重要的潮流计算方法，通过迭代计算逼近非线性方程组的解，可以得到电力系统中各节点的电压幅值和相角，用于分析电力系统的稳定性和安全性。
电力系统潮流计算完整程序及详细理论说明
电力系统潮流计算完整程序及详细理论说明——秦羽风在我刚开始学习潮流程序时,总是找不到一个正确的程序开始模仿学习。
后来经过多方努力，终于自己写出了一个结构清晰、完整的潮流程序。
此程序是一个通用的程序，只需要修改输入数据的子函数（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。
电力系统稳态分析--潮流计算
电力系统稳态分析--潮流计算(总36页)--本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--电力系统稳态分析摘要电力系统潮流计算是研究电力系统稳态运行情况的一种重要的分析计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各母线的电压，各元件中流过的功率，系统的功率损耗。
所以，电力系统潮流计算是进行电力系统故障计算，继电保护整定，安全分析的必要工具。
本文介绍了基于MATLAB软件的牛顿-拉夫逊法和P-Q分解法潮流计算的程序，该程序用于计算中小型电力网络的潮流。
在本文中，采用的是一个5节点的算例进行分析，并对仿真结果进行比较，算例的结果验证了程序的正确性和迭代法的有效性。
关键词：电力系统潮流计算；MATLAB；牛顿-拉夫逊法；P-Q分解法；目次1 绪论 ..................................................................................................... 错误!未定义书签。
背景及意义......................................................................................... 错误!未定义书签。
相关理论 ............................................................................................ 错误!未定义书签。
本文的主要工作 ................................................................................ 错误!未定义书签。
2 潮流计算的基本理论 ......................................................................... 错误!未定义书签。
电力系统潮流计算
(k ) f ( x ) (k ) x f ( x ( k ) )
迭代过程的收敛判据为 f ( x ( k ) ) 1
x ( k ) 2
或
牛顿—拉夫逊法实质上就是切线法，是一种逐步线性化的 方法。牛顿法不仅用于求解单变量方程，它也是求解多变 量非线性方程的有效方法。
有
(0) (0) (0) (0) f1 ( x1(0) x1(0) , x2 x2 , , xn xn )0 (0) (0) (0) (0) f 2 ( x1(0) x1(0) , x2 x2 , , xn xn )0
(0) (0) (0) (0) f n ( x1(0) x1(0) , x2 x2 , , xn xn )0
牛顿-拉夫逊法潮流计算
一、牛顿—拉夫逊法的基本原理 单变量非线性方程： x=x(0)+ Δx(0) 即 f(x=x(0)+ Δx(0) ) = 0 f(x)=0 (11—29) 解的近似值x(0)，它与真解的误差为Δx(0)
展成泰勒级数
f (x
(0)
x ) f ( x ) f ( x )x
f1 (0) xn )0 xn 0 f (0) 2 xn )0 xn 0
(0) f n ( x1(0) , x2 ,
写成矩阵形式：
f n f (0) x1(0) n x2 x1 0 x2 0 f1 x1 0 (0) (0) (0) f1 ( x1 , x2 , , xn ) f 2 (0) (0) (0) f 2 ( x1 , x2 , , xn ) x 1 0 (0) (0) (0) f ( x , x , , x n 1 2 n ) f n x1 0
5节点电力系统潮流计算matlab程序汇总
%1.形成节点导纳矩阵, yb55=6.250-18.750j;yb51=-5.000+15.000j;yb52=-1.250+3.750j;yb53=0.000-0.000j;yb54=0.000-0.000j; yb15=-5.000+15.000j;yb11=10.834-32.500j;yb12=-1.667+5.000j;yb13=-1.667+5.000j;yb14=-2.500+7.500j; yb25=-1.250+3.750j;yb21=-1.667+5.000j;yb22=12.917-38.750j;yb23=-10.000+30.000j;yb24=0.000-0.000j; yb35=0.000-0.000j;yb31=-1.667+5.000j;yb32=-10.000+30.000j;yb33=12.917-38.750j;yb34=-1.250+3.750j; yb45=0.000-0.000j;yb41=-2.500+7.500j;yb42=0.000-0.000j;yb43=-1.250+3.750j;yb44=3.750-11.250j; YB=[yb11 yb12 yb13 yb14 yb15; yb21 yb22 yb23 yb24 yb25 ;yb31 yb32 yb33 yb34 yb35; yb41 yb42 yb43 yb44 yb45 ;yb51 yb52 yb53 yb54 yb55]; disp('节点导纳矩阵YB=');disp(YB); %计算各节点功率的不平衡量设U=E+jF ；Y=G+Bj；E(1)=1.00;E(2)=1.00;E(3)=1.00;E(4)=1.00; F(1)=0;F(2)=0;F(3)=0;F(4)=0;G=real(YB);B=imag(YB); %设S=P+Bj； S(1)=0.20+0.20i;S(2)=-0.45-0.15i;S(3)=-0.40-0.05i;S(4)=-0.60-0.10i; P=real(S);Q=imag(S); k=0;precision=1; N1=4; while precision > 0.00001 E(5)=1.06;F(5)=0; for m=1:N1 for n=1:N1+1 %计算Pi,Qi,设Pi=Pt;Qi=Qt Pt(n)=(E(m)*(G(m,n)*E(n)-B(m,n)*F(n))+F(m)*(G(m,n)*F(n)+B(m,n)*E(n))); Qt(n)=(F(m)*(G(m,n)*E(n)-B(m,n)*F(n))-E(m)*(G(m,n)*F(n)+B(m,n)*E(n))); end %设P,Q的改变量为dP,dQ dP(m)=P(m)-sum(Pt); dQ(m)=Q(m)-sum(Qt); end for m=1:N1 for n=1:N1+1 %计算Hij Nij Jij Lij H(m,n)=-B(m,n)*E(m)+G(m,n)*F(m);N(m,n)=G(m,n)*E(m)+B(m,n)*F(m); J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m);L(m,n)=G(m,n)*F(m)-B(m,n)*E(m); end end for m=1:N1 for n=1:N1+1Bi(n)=G(m,n)*F(n)+B(m,n)*E(n); Ai(n)=G(m,n)*E(n)-B(m,n)*F(n); end %sum(Ai)，sum(Bi)用于实现公式中的sigerma从j到n的求和; H(m,m)=sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m))+2*G(m,m)*F(m); N(m,m)=sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m))+2*G(m,m)*E(m); J(m,m)=-2*B(m,m)*F(m)+sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m)); L(m,m)=-2*B(m,m)*E(m)-(sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m))); end %设雅可比矩阵为JJ，以下语句用来实现雅可比矩阵中对角线上元素H N J L 的排列 for m=1:N1 JJ(2*m-1,2*m-1)=H(m,m); JJ(2*m-1,2*m)=N(m,m); JJ(2*m,2*m-1)=J(m,m); JJ(2*m,2*m)=L(m,m); end %以下语句用于实现雅可比矩阵非对角线上元素的排列 for m=1:N1 for n=1:N1 if m==n else H(m,n)=-B(m,n)*E(m)+G(m,n)*F(m); N(m,n)=G(m,n)*E(m)+B(m,n)*F(m); J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m); L(m,n)=G(m,n)*F(m)-B(m,n)*E(m); JJ(2*m-1,2*n-1)=H(m,n); JJ(2*m-1,2*n)=N(m,n); JJ(2*m,2*n-1)=J(m,n); JJ(2*m,2*n)=L(m,n); end end end %设由P,Q的改变量组成的8×1矩阵为PQ，由E,F的改变量组成的8×1矩阵为dU for m=1:N1 PQ(2*m-1)=dP(m); PQ(2*m)=dQ(m); end dU=inv(JJ)*PQ';precision=max(abs(dU)); for n=1:N1 F(n)=F(n)+dU(2*n-1); E(n)=E(n)+dU(2*n); end for n=1:N1+1 U(n)=E(n)+(F(n))*j; end k=k+1; k-1, dU=dU',PQ,U end。
5节点电力系统牛顿-拉夫逊法潮流计算
（二 〇 一 四 年 十 二 月课 程 论 文 学校代码： 10128 学 号: 20141100304题 目：五节点系统计算机潮流计算编程 学生姓名：张佳羽学 院：电力学院系 别：电力系专 业：电力系统及其自动化指导教师：郭力萍程序设计% 本程序的功能是用牛顿拉夫逊法进行潮流计算n=input('请输入节点数:n=’);nl=input(’请输入支路数：nl='）；isb=input('请输入平衡母线节点号:isb=’);pr=input('请输入误差精度:pr=’）；B1=input（’请输入由各支路参数形成的矩阵：B1=')；B2=input(’请输入各节点参数形成的矩阵：B2=’）；X=input（'请输入由节点号及其对地阻抗形成的矩阵:X=’）；Y=zeros(n）；e=zeros(1，n);f=zeros（1，n);V=zeros（1，n）；O=zeros(1，n）；S1=zeros(nl);for i=1:nlif B1(i，6)==0p=B1（i，1);q=B1(i,2)；else p=B1（i,2);q=B1(i，1)；endY（p，q)=Y(p，q）—1。
/(B1（i,3）*B1（i,5））;Y(q,p）=Y（p，q）；Y(q，q）=Y(q，q）+1。
/(B1(i,3)＊B1（i,5)^2）+B1（i,4)./2；Y(p,p)=Y（p，p)+1。
/B1（i,3）+B1(i，4)./2;end%求导纳矩阵disp('导纳矩阵Y=')；disp(Y）;G=real(Y）；B=imag(Y）;for i=1：ne(i）=real(B2（i,3))；f(i)=imag（B2（i，3));V(i)=B2(i,4）;endfor i=1:nS（i）=B2（i，1)—B2（i,2）；B(i,i)=B(i，i)+B2（i,5）；endP=real(S);Q=imag(S)；ICT1=0；IT2=1;N0=2*n;N=N0+1;a=0；while IT2~=0IT2=0；a=a+1;for i=1：nif i~=isbC（i）=0；D（i）=0；for j1=1：nC（i）= C(i）+G（i，j1）*e（j1)—B(i，j1)＊f(j1）;D(i)= D（i）+G(i,j1)*f(j1）+B(i,j1)*e(j1);endP1=C（i）*e(i)+f（i)＊D(i);Q1=f（i)＊C(i)—D(i)*e(i);V2=e(i）^2+f(i）^2;if B2（i,6）~=3DP=P（i)-P1;DQ=Q(i)-Q1;for j1=1：nif j1~=isb＆j1~=iX1=-G(i，j1）*e(i）—B（i,j1）*f(i）;X2=B(i,j1）*e(i)—G(i，j1）＊f(i);X3=X2;X4=-X1；p=2＊i-1；q=2*j1—1；J（p，q）=X3;J（p,N)=DQ;m=p+1;J(m,q)=X1；J（m,N）=DP;q=q+1;J(p,q)=X4；J（m，q）=X2;elseif j1==i&j1～=isbX1=—C（i）-G(i，i）*e（i)—B（i,i）*f(i);X2=-D(i)+B(i,i)*e(i)-G（i，i）＊f（i）；X3=D（i)+B（i，i)*e(i）-G（i,i）*f(i）；X4=-C(i)+G（i,i）*e（i）+B(i，i）*f（i）；p=2*i-1;q=2＊j1—1；J(p，q）=X3;J（p,N)=DQ;m=p+1;J(m，q)=X1;J（m，N)=DP；q=q+1;J(p，q)=X4;J(m，q）=X2;endendelseDP=P(i）—P1;DV=V(i）^2—V2；for j1=1：nif j1～=isb＆j1～=iX1=—G(i，j1)*e(i）-B(i，j1)＊f（i）；X2=B（i,j1)*e（i)—G（i,j1）*f（i）；X5=0;X6=0；p=2*i-1;q=2＊j1-1；J(p,q)=X5；J(p，N)=DV；m=p+1；J(m，q)=X1；J(m，N)=DP；q=q+1；J(p,q）=X6;J（m，q）=X2;elseif j1==i＆j1～=isbX1=—C(i）—G（i，i）＊e(i）—B（i，i）*f(i）；X2=-D(i）+B(i，i)＊e（i)-G(i，i)＊f（i)；X5=-2＊e(i）；X6=—2*f（i);p=2*i-1；q=2*j1—1；J（p，q)=X5;J（p,N）=DV;m=p+1；J(m，q)=X1；J(m，N)=DP；q=q+1；J(p,q）=X6；J(m，q）=X2;endendendendend％求雅可比矩阵for k=3：N0k1=k+1;N1=N；for k2=k1：N1J(k,k2)=J(k,k2)./J(k，k）；endJ(k,k)=1;if k~=3;k4=k—1；for k3=3:k4for k2=k1：N1J（k3,k2）= J(k3，k2)—J(k3,k)＊J(k，k2）；endJ(k3，k)=0；endif k==N0，break；endfor k3=k1:N0for k2=k1:N1J(k3，k2)=J(k3,k2)-J(k3，k）*J(k,k2);endJ（k3,k）=0；endelsefor k3=k1：N0for k2=k1:N1J（k3,k2）= J（k3,k2）—J（k3,k)*J（k，k2)；endJ（k3，k）=0;endendendfor k=3：2:N0—1L=(k+1)。
潮流计算步骤
潮流计算步骤
潮流计算是电力系统分析中的一种基本计算方法，用于确定电网中的电压分布和功率流动情况。
以下是潮流计算的基本步骤：
1、输入原始数据和信息：包括电网的结构信息、设备参数、负荷和电源的分布及大小等。
2、建立数学模型：根据电路理论和电力系统网络模型，建立描述电力系统中电压、电流和功率关系的数学模型。
3、形成节点导纳矩阵：根据电网结构，形成节点导纳矩阵，用于描述系统中各节点之间的电气联系。
4、确定待求状态变量初值：根据实际情况，为待求的状态变量（如节点电压）设定初值。
5、迭代求解：使用迭代法对数学模型进行求解，逐步更新状态变量的值，直到满足收敛条件为止。
6、计算节点电压：根据迭代求解的结果，计算出各节点的电压值。
7、计算功率分布：根据节点电压和网络参数，计算出各支路的功率流动情况。
8、结果分析：对计算结果进行整理和分析，评估电网的运行状态，为进一步优化和调整提供依据。
需要注意的是，潮流计算的具体步骤可能会因不同的计算方法和电力系统分析软件而有所差异。
在实际应用中，需要根据具体的软件
和要求进行操作。
电力系统潮流计算程序(详细)
潮流计算，电力学名词，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。
潮流计算是根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。
通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。
待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。
潮流计算程序，相关的原始数据数据数据输入格式如下：%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。
牛拉法潮流计算程序（附3机9节点结果对比）
摘要电力系统潮流计算是研究电力系统稳态运行的一种重要方法,它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态,包括各母线的电压、线路的功率分布以及功率损耗等等。
潮流计算主要用于电网规划和静态安全分析，它可为扩建电力网络，以达到规划周期内所需要的输电能力提供依据；也可以对预想事故进行模拟和分析，校核预想事故下的电力系统安全性。
本文简单介绍了牛顿-拉夫逊潮流计算的原理、模型与算法，然后用具体的实例，利用MATLAB对牛顿-拉夫逊法的算法进行了验证。
关键词：电力系统潮流计算牛顿-拉夫逊法 MATLAB一、牛拉法的数学模型对一个N 节点的电力网路，列写节点电压方程，即I =Y V(1.1)式中，I 为节点注入电流列相量，Y 为节点导纳矩阵，V 为节点电压列相量。
由于异地测量的两个电流缺少时间同步信息，以注入功率替换注入电流作为已知量。
即***1+niij j ij j i i i Y V V I V Q P ∙∙===∑(1.2)其中,Y ij =G ij +j B ij ，带入上式，得到有功功率和无功功率方程 P i =V i V j G ij cos θij +B ij sin θij n j =1 (1.3)Q i =V i V j G ij sin θij −B ij cos θij n j =1(1.4)大部分情况下，已知PQ ，求解V θ。
考虑到电网的功率平衡，至少选择一台发电机来平衡全网有功功率，即至少有一个平衡节点，常选择调频或出线较多的发电机作为平衡节点。
具有无功补偿的母线能保持电压幅值恒定，这类节点可作为PV 节点。
潮流计算中节点分类总结如下：已知电力系统有m 个PQ 节点，r 个PV 节点和1个平衡节点，则可以提取m+r 个有功功率方程和m 个无功功率方程，从而求解出m+r 个θ和m 个V ，其余节点的有功和无功可通过式（1.3）、（1.4）求得，这样就完成了潮流计算。
两机五节点网络潮流计算
两机五节点网络潮流计算潮流计算是一种基本电力系统分析方法，用于确定电力系统中各个节点的电压幅值和相位角，以及线路中的功率流动情况。
潮流计算的目的是为了评估电力系统的稳定性并优化系统操作。
在两机五节点网络中，可以假设两个发电机节点为节点1和2，其他五个节点为节点3、4、5、6和7、根据电力系统的潮流计算原理，可以按照下面的步骤进行计算：步骤1：定义网络拓扑和参数首先，需要定义网络的拓扑结构和电气参数。
拓扑结构指的是节点之间的连接关系，参数包括电压幅值、相位角、阻抗等。
步骤2：初始化网络状态初始化电压和功率的初值。
可以假设所有节点的电压初始值为1、假设发电机节点1的有功功率为1，节点2的有功功率为0。
步骤3：计算节点注入功率根据节点电压和电流的关系，可以计算出每个节点的注入功率。
对于发电机节点，注入功率可以通过给定的有功功率来计算。
对于非发电机节点，注入功率可以通过节点电压、导纳矩阵和节点电流来计算。
步骤4：更新节点电压根据节点注入功率和导纳矩阵，可以计算出每个节点的电流注入。
根据电流和节点电压的关系，可以更新节点的电压。
步骤5：迭代计算重复步骤3和步骤4，直到达到收敛条件。
收敛条件通常是指节点电压的变化小于一个给定的阈值，如0.001步骤6：计算线路功率在达到收敛状态后，可以根据节点电压和导纳矩阵计算出线路的功率。
潮流计算的结果是电力系统中各个节点的电压幅值和相位角，以及线路中的功率流动情况。
这些结果对于电力系统的运行和规划非常重要，可以评估系统的稳定性，并为系统操作和优化提供依据。
总之，两机五节点网络潮流计算是电力系统潮流计算的一种常见场景。
通过定义网络拓扑和参数，初始化网络状态，计算节点注入功率，更新节点电压，迭代计算，最后计算线路功率，可以得到电力系统中各个节点的电压和功率的结果。
这些结果对于电力系统的运行和规划非常重要。
潮流计算程序及计算结果
附表1：计算机计算潮流程序：%本程序的功能是用牛顿——拉夫逊法进行潮流计算% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳% 5、支路的变比；6、支路首端处于K侧为1，1侧为0% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量% 6、节点分类标号clear;n=13;%input('请输入节点数：n=');nl=13;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[];%input('请输入由支路参数形成的矩阵：B1=');B2=[];%input('请输入各节点参数形成的矩阵：B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl); %－－－－－－－修改部分－－－－－－－－－－－－ym=0;SB=100;UB=220;%ym=input('您输入的参数是标么值？（若不是则输入一个不为零的数值）'); if ym~=0%SB=input('请输入功率基准值:SB=');%UB=input('请输入电压基准值:UB=');YB=SB./UB./UB;BB1=B1;BB2=B2;for i=1:nlB1(i,3)=B1(i,3)*YB;B1(i,4)=B1(i,4)./YB;enddisp('B1矩阵B1=');disp(B1)for i=1:nB2(i,1)=B2(i,1)./SB;B2(i,2)=B2(i,2)./SB;B2(i,3)=B2(i,3)./UB;B2(i,4)=B2(i,4)./UB;B2(i,5)=B2(i,5)./SB;enddisp('B2矩阵B2=');disp(B2)end% % %---------------------------------------------------for i=1:nl %支路数if B1(i,6)==0 %左节点处于低压侧p=B1(i,1);q=B1(i,2);elsep=B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非对角元Y(q,p)=Y(p,q);Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %对角元K侧Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元1侧end%求导纳矩阵disp('导纳矩阵Y=');disp(Y)%----------------------------------------------------------G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部for i=1:n %给定各节点初始电压的实部和虚部e(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4); %PV节点电压给定模值endfor i=1:n %给定各节点注入功率S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SLB(i,i)=B(i,i)+B2(i,5); %i节点无功补偿量end%=========================================================== ========P=real(S);Q=imag(S);ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;while IT2~=0IT2=0;a=a+1;for i=1:nif i~=isb %非平衡节点C(i)=0;D(i)=0;for j1=1:nC(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)endP1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求P',Q'V2=e(i)^2+f(i)^2; %电压模平方%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素=========if B2(i,6)~=3 %非PV节点DP=P(i)-P1; %节点有功功率差DQ=Q(i)-Q1; %节点无功功率差%=============== 以上为除平衡节点外其它节点的功率计算=================%================= 求取Jacobi矩阵===================for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de=-dQ/dfX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/df=dQ/deX3=X2; % X2=dp/df X3=dQ/deX4=-X1; % X1=dP/de X4=dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;elseif j1==i&j1~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dQ/deX4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Qm=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△PJ(m,q)=X2;endendelse%=============== 下面是针对PV节点来求取Jacobi矩阵的元素===========DP=P(i)-P1; % PV节点有功误差DV=V(i)^2-V2; % PV节点电压误差for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/deX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/dfX5=0;X6=0;p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;elseif j1==i&j1~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX5=-2*e(i);X6=-2*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;endendendendend%========= 以上为求雅可比矩阵的各个元素===================== for k=3:N0 % N0=2*n （从第三行开始，第一、二行是平衡节点）k1=k+1;N1=N; % N=N0+1 即N=2*n+1扩展列△P、△Qfor k2=k1:N1 % 扩展列△P、△QJ(k,k2)=J(k,k2)./J(k,k); % 非对角元规格化endJ(k,k)=1; % 对角元规格化if k~=3 % 不是第三行%======================================================== ====k4=k-1;for k3=3:k4 % 用k3行从第三行开始到当前行前的k4行消去for k2=k1:N1 % k3行后各行下三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endif k==N0break;end%==========================================for k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endelsefor k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endendend%====上面是用线性变换方式将Jacobi矩阵化成单位矩阵=====for k=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N); %修改节点电压实部k1=k+1;f(L)=f(L)-J(k1,N); %修改节点电压虚部end%------修改节点电压-----------for k=3:N0DET=abs(J(k,N));if DET>=pr %电压偏差量是否满足要求IT2=IT2+1; %不满足要求的节点数加1endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解"w=-J*V"disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:nV(k)=sqrt(e(k)^2+f(k)^2);sida(k)=atan(f(k)./e(k))*180./pi;E(k)=e(k)+f(k)*j;end%=============== 计算各输出量=========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);EE=E*UB;disp(EE);disp('-----------------------------------------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);VV=V*UB;disp(VV);disp('-----------------------------------------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);for p=1:nC(p)=0;for q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));endS(p)=E(p)*C(p);enddisp('各节点的功率S为(节点号从小到大排列)：');disp(S);disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');SS=S*SB;disp(SS);disp('-----------------------------------------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./( B1(i,3)*B1(i,5))));Siz(i)=Si(p,q);elseSi(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))-conj(E(q)))*conj(1./( B1(i,3)*B1(i,5))));Siz(i)=Si(p,q);enddisp(Si(p,q));SSi(p,q)=Si(p,q)*SB;ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);%disp(SSi(p,q));disp('-----------------------------------------------------');enddisp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./( B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);elseSj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))-conj(E(p)))*conj(1./( B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);enddisp(Sj(q,p));SSj(q,p)=Sj(q,p)*SB;ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);%disp(SSj(q,p));disp('-----------------------------------------------------');enddisp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);disp(DS(i));DDS(i)=DS(i)*SB;ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);%disp(DDS(i));disp('-----------------------------------------------------');endfigure(1);subplot(2,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(2,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;subplot(2,2,3);bar(real(S));ylabel('节点注入有功');grid on;subplot(2,2,4);bar(Siz);ylabel('支路首端无功');grid on;1.冬季最大运行方式潮流计算结果:计算机运行的B1,B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.88+0.545*i 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.88+0.545*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.2+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表：节点号电压值相角值支路标号首端功率末端功率支路功率损耗1 242.0000 0 1-2 148.9391+3.327637i -143-27.45476i 5.93913-24.1271i2 231.0000 -3.0341 1-11 89.1956+37.193i -88.19803-61.4687i 0.997519-24.2757i3 227.4154 -4.4383 11-3 88.19803+61.4687i -88-54.5i 0.19803+6.9687i4 226.3304 -4.9004 1-12 77.8364+36.4186i -77.2404-55.7062i 0.596039-19.2876i5 229.8318 -3.9547 12-4 77.2404+55.7062i -77-47.72i 0.24036+7.9862i6 217.7226 -8.2720 1-5 63.5438+12.2204i -61.8773-32.0321i 1.66656-19.8116i7 234.9326 -3.1047 5-6 77.2597+56.3501i -77-47.72i 0.25974+8.6301i8 226.0991 -6.2429 5-7 15.3825-24.3181i 15.5354+0.274108i 0.152961-24.044i9 231.0000 0.5851 7-8 88.20085+61.55009i -88-54.5i 0.20085+7.0501i10 231.0000 3.4950 1-7 40.806-6.05737i -40.0647-22.0555i 0.741264-28.1128i11 236.1931 -1.3348 7-13 -63.6716-39.7687i 64.0965+17.5832i 0.424934-22.1855i12 237.9346 -0.8892 13-9 -64.0965-17.5832i 64.2+21.0187i 0.10348+3.4355i13 238.1868 -2.2028 1-10 79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果图形：2.冬季最小运行方式潮流计算结果：计算机运行的B1B2矩阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.616+0.3817*i 1 0 0 20 0.539+0.3817*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.1+0.06197*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]电压调整措施：变电所1、4变压器变比：+2.5% 水电厂变压器变比：-2.5%5 234.2241 -2.0051 (12--4) 54.0234+42.2706i -53.9-38.17i 0.12342+4.1006i6 226.2159 -4.8577 (1--5) 34.1669+4.22856 -33.6427-28.276i 0.524165-24.0474i7 235.1128 -0.4737 (5--6) 54.0179+37.3169i -53.9-33.4i 0.11789+3.9169i8 223.9941 -2.4603 (5--7) -20.3752-9.04094i 20.5371-15.4558i 0.161947-24.4968i9 231.0000 3.4429 (1--7) 11.0013+1.45186i -10.8258-31.4513i 0.175442-29.9995i10 231.0000 3.9321 (7--8) 53.9768+36.0957i -53.9-33.4i 0.076797+2.6957i11 238.4187 -0.9561 (7--13) -63.6881+10.8115i 64.0874-32.7763i 0.39932-21.9648i12 239.1742 -0.6185 (13--9) -64.0874+32.7763i 64.2-29.0386i 0.11258+3.7377i13 234.9468 0.6942 (1--10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机运行结果的图形：3.夏季最大运行方式计算机计算结果：计算机运行B1B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.776+0.481*i 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.776+0.481*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表：10 231.0000 3.4950 (7,8) 77.7585+53.6634i -77.6-48.1i0.1585+5.5634i11 237.0938 -1.1858 (7,13) -113.4984+9.476406i 114.6926-28.38885i 1.19422-18.9124i12 239.1742 -0.6185 (13,9) -114.6926+28.38885 115-18.18369i0.307384+10.2052i13 233.3912 3.2303 (1,10) -79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果如图：4.夏季最小运行方式：计算机运行B1B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.543+0.336*i 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.543+0.336*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]10 231.0000 3.9321 (7,8) 54.3735+36.2509i -54.3-33.67i 0.073528+2.5809i11 239.0099 -0.8518 (7,13) -113.4428+24.39707i 114.6754-43.79301i 1.23255-19.3959i12 239.8726 -0.5689 (13,9) -114.6754+43.79301i 115-33.01613i 0.324605+10.7769i13 235.9565 4.4510 (1,10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机计算结果如图：5.夏季故障运行状态:调压及无功补偿措施如下:变电所3的变压器变比为-2.5%,无功补偿容量为20Mvar。
电力系统潮流计算C程序精编版
电力系统潮流计算注：这是一个基于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]);}}。
电力系统潮流及短路电流计算程序
电力系统潮流及短路电流计算程序以下是一个简单的电力系统潮流计算程序的框架：1.输入数据准备阶段：-输入潮流计算的电力系统拓扑结构，包括各节点之间的连接关系、导线电阻、电抗等信息。
-输入电力系统的负荷信息，包括负荷节点、负荷大小、负荷类型等。
-输入电力系统的发电机信息，包括发电机节点、发电机类型、发电机容量等。
2.潮流计算阶段：-初始条件设置：给定电力系统中各节点的初始电压、相角等信息。
-节点功率方程求解：根据电力系统的拓扑结构和发电机、负荷信息，建立节点功率方程。
-潮流计算迭代：使用牛顿-拉夫逊法等迭代算法求解节点功率方程，得到各节点的电压、相角等参数。
3.潮流计算结果输出阶段：-输出各节点的电压、相角、有功功率、无功功率等参数。
-输出各支路的电流、功率损耗等参数。
-输出系统的功率平衡情况。
4.短路电流计算阶段：-输入短路电流计算的电力系统拓扑结构。
-输入短路电流计算的负荷信息。
-输入短路电流计算的电源信息。
-使用KVL（电压法）或KCL（电流法）等方法计算各节点短路电流。
5.短路电流计算结果输出阶段：-输出各节点的短路电流大小。
-输出各支路的短路电流大小。
以上只是一个电力系统潮流及短路电流计算程序的大致流程框架，具体实现细节和算法选择还需要根据具体情况进行进一步的设计和开发。
在实际应用中，还需要考虑各种特殊情况和计算优化方法，以提高计算速度和准确性。
总之，电力系统潮流及短路电流计算程序是电力工程师在设计和运行电力系统中不可或缺的工具，它能够帮助工程师快速了解系统的运行状态和电流分布情况，以便进行系统优化和安全评估。
电力系统潮流计算实验报告
１1． 手算过程已知：节点1：PQ 节点， s(1)= -0.5000-j0.3500 节点2：PV 节点， p(2)=0.4000 v(2)=1.0500 节点3：平衡节点，U(3)=1.0000∠0.0000 网络的连接图：0.0500+j0.2000 1 0.0500+j0.2000231）计算节点导纳矩阵由2000.00500.012j Z 71.418.112j y ;2000.00500.013j Z71.418.113j y ;导纳矩阵中的各元素：42.936.271.418.171.418.1131211j j j y y Y ;71.418.11212j y Y ; 71.418.11313j y Y; 21Y 71.418.11212j y Y ; 71.418.12122j y Y;002323j y Y;31Y 71.418.11313j y Y; 32Y 002323j y Y;71.418.13133j y Y;形成导纳矩阵BY ：71.418.10071.418.10071.418.171.418.171.418.171.418.142.936.2j j j j j j j j j Y B2）计算各PQ、PV 节点功率的不平衡量，及PV 节点电压的不平衡量：取:000.0000.1)0(1)0(1)0(1j jf e U000.0000.1)0(2)0(2)0(2j jf e U节点3是平衡节点，保持000.0000.1333j jf e U为定值。
nj j jij jij ijij jij i ieB fG f fB eG e P1)0()0()0()0()0()0()0(；２nj j jij jij ijij jij i ie B fG e f B eG f Q 1)0()0()0()0()0()0()0(；);(2)0(2)0(2)0(iiif e U)0.142.90.036.2(0.0)0.042.90.136.2(0.1)0(1P)0.171.40.018.1(0.0)0.071.40.118.1(0.1 )0.171.40.018.1(0.0)0.071.40.118.1(0.1 0.0 ；)0.142.90.036.2(0.1)0.042.90.136.2(0.0)0(1Q)0.171.40.018.1(0.1)0.071.40.118.1(0.0 )0.171.40.018.1(0.1)0.071.40.118.1(0.0 0.0 ；)0.171.40.018.1(0.0)0.071.40.118.1(0.1)0(2P)0.171.40.018.1(0.0)0.071.40.118.1(0.1 )0.00.00.00.0(0.0)0.10.00.10.0(0.1 0.0 ；101)(222)0(22)0(22)0(2f e U；于是：;)0()0(iiiP P P ;)0()0(iiiQQ Q);(2)0(2)0(22)0(iiiif e UU5.00.05.0)0(11)0(1P P P ;35.00.035.0)0(11)0(1QQ Q;4.00.04.0)0(22)0(2P P P ;1025.0)01(05.1)(2222)0(22)0(2222)0(2f e UU3）计算雅可比矩阵中各元素雅可比矩阵的各个元素分别为：３ji ij ji ij j i ij j i ij ji ij j i ij e U S f U R e Q L f Q J e P N f P H 22;;; 又： nj j jij jij i jij jij i ieB fG f fB eG e P1)0()0()0()0()0()0()0(； nj j jij jij ijij jij iieB fG e fB eG f Q 1)0()0()0()0()0()0()0(；);(2)0(2)0(2)0(iiif e U)0(1P )0(111)0(111)0(1)0(111)0(111)0(1e Bf G f f B e G e)0(212)0(212)0(1)0(212)0(212)0(1e B fG f f B e G e313313)0(1313313)0(1e Bf G f f B e G e ；)()()0(111)0(111)0(1)0(111)0(111)0(1)0(1e Bf Ge f B e G f Q)()()0(212)0(212)0(1)0(212)0(212)0(1e Bf G e f B e G f)()(313313)0(1313313)0(1e Bf G e f B e G f；)0(2P )0(121)0(121)0(2)0(121)0(121)0(2e Bf G f f B e G e)0(222)0(222)0(2)0(222)0(222)0(2eB fG f fBeG e323323)0(2323323)0(2e Bf G f f B e G e ；)(2)0(22)0(22)0(2f e U42.90.171.40.171.4313)0(212)0(1)0(1)0(11e B e Bf P H ； 36.20.118.10.118.10.136.222313)0(212)0(111)0(1)0(1)0(11 e G e G e G e P N 36.20.118.10.118.1313)0(212)0(1)0(1)0(11 e G e G f Q J４42.90.171.40.171.40.142.922313)0(212)0(111)0(1)0(1)0(11 e B e B e B e Q L 71.40.171.4)0(112)0(2)0(1)0(12 e B f P H ； 18.10.118.1)0(112)0(2)0(1)0(12 e G e P N ； 18.10.118.1)0(112)0(2)0(1)0(12 e G f Q J ；71.40.171.4)0(112)0(2)0(1)0(12 e B e Q L ； 71.40.171.4)0(221)0(1)0(2)0(21 e B f P H ； 11.40.111.4)0(221)0(1)0(2)0(21 e G e P N ； 0)0(12)0(2)0(21 f U R ； 0)0(12)0(2)0(21 e U S ； 71.40.10.00.171.4323)0(121)0(2)0(2)0(22 e B e B f P H ； 18.10.10.00.118.10.118.122323)0(121)0(222)0(2)0(2)0(22 e G e G e G e P N ；02)0(2)0(22)0(2)0(22 f f U R ； 0.20.122)0(2)0(22)0(2)0(22 e e U S ； 得到K=0时的雅可比矩阵：0.200018.171.418.171.471.418.142.936.218.171.436.242.9)0(J4）建立修正方程组：５)0(2)0(2)0(1)0(10.200011.4959.1011.4959.10959.1011.4918.2122.811.4959.1022.8918.210975.04.035.08.0e f e f 解得：04875.001828.00504.00176.0)0(2)0(2)0(1)0(1e f e f 因为 )0()0()1(iiie e e ； )0()0()1(iiif f f ；所以 9782.00218.00.1)0(1)0(1)1(1e e e ； 0158.00158.00)0(1)0(1)1(1f f f ；05125.105125.00.1)0(2)0(2)1(2e e e ；05085.005085.00)0(2)0(2)1(2f f f ；5）运用各节点电压的新值进行下一次迭代：即取: 0158.09782.0)1(1)1(1)1(1j jf e U05085.005125.1)1(2)1(2)1(2j jf e U节点3时平衡节点，保持000.0000.1333j jf e U为定值。
ieee5节点标准测试系统潮流计算与仿真原理
ieee5节点标准测试系统潮流计算与仿真原理IEEE 5节点标准测试系统是电力系统研究中使用最广泛的标准测试系统之一。
该系统由5个节点组成，包括1个发电机节点、3个负荷节点和1个传输线节点。
这一简化的系统用于进行各种电力系统潮流计算和仿真研究，能够帮助我们理解电力系统的性能和行为。
在电力系统的潮流计算和仿真中，我们需要使用各种数学方法和计算技术来求解系统中的电压、电流和功率等参数。
这些参数对于电力系统的可靠运行和稳定性至关重要。
潮流计算是电力系统分析中最基本的问题之一，它通过求解节点电压和功率的方程，来确定电力系统中各个元件之间的功率分配情况。
潮流计算的结果可以用于评估系统的稳态损耗、潮流分布、功率负荷和发电机容量等信息。
IEEE 5节点标准测试系统的潮流计算和仿真可以通过直接计算法和迭代法两种方法进行。
直接计算法是根据系统拓扑结构和各个元件的参数，直接推导潮流方程并求解。
这种方法适用于计算较小规模的电力系统，但对于大规模系统计算量较大。
迭代法则是通过不断迭代调整节点电压和功率的方法来求解潮流方程。
迭代法需要通过选择合适的迭代算法和初始猜测值来保证收敛性，但适用于计算较大规模的系统。
IEEE 5节点标准测试系统的潮流计算和仿真需要构建系统的数学模型。
该系统中的节点可以表示为复数形式的变量，节点之间的链接可以表示为复数形式的导纳。
潮流计算中的潮流方程可以表示为以下形式：I = Y * V其中，I为节点电流，Y为节点导纳矩阵，V为节点电压。
通过将节点电压和功率作为未知数，将导纳矩阵作为已知量，可以将潮流方程转化为线性方程组，通过求解线性方程组即可得到节点电压和功率的值。
IEEE 5节点标准测试系统的潮流计算和仿真还可以用于研究系统的稳定性和可靠性。
通过改变系统中的负荷、发电机容量和传输线参数等，可以模拟不同的运行情况，并分析系统在不同负荷和故障条件下的稳定性和可靠性。
通过潮流计算和仿真，我们可以验证电力系统的稳态约束和安全边界，并针对性地优化系统的配置和控制策略。
两机五节点网络潮流计算—牛拉法
两机五节点网络潮流计算—牛拉法牛拉法（Gauss-Seidel Method）是一种常用的迭代方法，用于解决电力系统的潮流计算问题。
在电力系统中，潮流计算是一项重要的工作，用于求解网络中各节点的电压和功率大小。
牛拉法是一种有效的求解方法，适用于小型电力系统，其基本思想是通过迭代来逼近最优解。
潮流计算问题可以抽象成求解非线性方程组的问题，即求解节点电压复数值的方程组。
具体来说，我们需要求解以下方程组：P_i = V_i * ( G_ii * cosθ_i + ∑(G_ij * cos(θ_i - θ_j)) - B_ii * sinθ_i - ∑(B_ij * sin(θ_i - θ_j)))Q_i = V_i * ( G_ii * sinθ_i + ∑(G_ij * sin(θ_i - θ_j)) + B_ii * cosθ_i + ∑(B_ij * cos(θ_i - θ_j)))其中，P_i和Q_i分别表示第i个节点的有功功率和无功功率，V_i表示第i个节点的电压幅值，θ_i表示第i个节点的电压相角，G_ij和B_ij分别表示节点i和节点j之间的导纳和电纳。
牛拉法的基本思路是通过迭代，逐步逼近节点电压的最优解。
假设我们需要求解的是一个两机五节点网络。
首先，我们可以随机初始化每个节点的电压幅值和相角值（也可以根据经验给定初始值）。
然后，根据上述方程组，计算每个节点的有功功率和无功功率。
接下来，我们采用牛拉法的迭代步骤来逼近节点电压的最优解。
具体步骤如下：1.选择一个初始节点（可以是任意节点），将其电压相角θ_i固定为0。
2.通过方程组计算该节点的电压幅值V_i。
3.将计算得到的电压幅值V_i和电压相角θ_i作为该节点的新的电压值。
4.对于其他节点，计算它们的电压相角θ_i和电压幅值V_i，并将其更新为新的电压值。
5.重复2-4步骤，直到收敛或满足收敛条件。
在每次迭代过程中，我们可以根据收敛准则来判断是否达到收敛，通常是通过计算两次迭代之间电压的变化量来判断。
电力系统之五节点潮流计算
辽宁工程技术大学电力系统分析综合训练一设计题目 5节点电力网络潮流计算指导教师院（系、部）电气与控制工程学院专业班级智能电网信息工程学号 **********姓名日期 2016/05/04智能电网系综合训练标准评分模板说明：1. 线路1为双回线路，双回线路参数完全相同，其余线路为单回线路；2. 变压器T2为两台并联，并联运行的变压器参数相同，其他变压器为单台运行。
目录1本次综合训练目的： (1)2 Power World软件简介： (1)3单线图： (3)4手工计算导纳矩阵，并与软件计算结果比较： (5)5单步运行牛-拉法潮流计算，给出前3步计算结果： (6)6 在保证每条线路和变压器不过载的情况下，确定节点3上发电机的允许出力范围： (9)7 在节点2上添加200Mvar并联电容器组，观察节点2电压变化和线路损耗： (10)8静态安全分析：对“单个线路”和“单个变压器”进行N-1校验分析。
(11)9结论： (12)1本次综合训练目的：通过对多节点电网的潮流计算，巩固和运用前面所学到的潮流计算基础理论知识，掌握电力系统潮流计算机计算的一般原则和方法，掌握潮流计算软件的使用方法，培养学生分析问题和解决问题的能力。
2 Power World 软件简介：Power World Simulator 是一款电力系统仿真软件。
核心是：全面、强大的潮流计算程序，可有效计算100000个节点的电力网络。
工具栏上的两种模式：编辑模式（建立、修改实例和单线图）。
运行模式（对实例进行各种仿真计算）。
在编辑模式的绘图功能区中可以建立单线图，并设置修改单线图中各元件的参数。
Power World 的功能有基本的潮流计算、短路计算、静态安全分析、灵敏度分析、分布因子计算、分时段仿真、最优潮流、计及安全约束的最优潮流，暂态稳定计算。
潮流计算方法简介：电力系统潮流计算时研究电力系统稳定运行情况的一种基本电气计算。
他的任务是根据给定的运行条件和网络结构确定整个系统的运行状态，如母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
% A为支路参数形成的矩阵;
% (1)为该支路首端节点号;
% (2)为该支路首端节点号;
% (3)为支路阻抗;
%(4)为支路变压器变比;
A=[2,3,0.145+0.581i,0;
1,5,0.104+0.581i,0;
3,5,0.031+0.248i,0;
3,4,0.082+0.427i,0.95;
input('系统数据=');
Y=zeros(n);
%把发电机电抗加入自阻抗.....................................................
for i=1:n
if X(i,2)~=0;
p=X(i,1);
Y(p,p)=1./X(i,2);
end
end
%算节点阻抗矩阵............................................................
[l m]=size(A);
for i=1:l;
p=A(i,1);
q=A(i,2);
Y(p,q)=-1./A(i,3);
f=3;
%故障分量网络节点电压方程，阻抗形式..........................................
Zf=zeros(n,1);
for i=1:n
Zf(i)=Z(i,f);
end %把短路支路取一列阻抗...................................................
disp(IG1);
IG2=0.5*(-dU(2))/X(2,2);
disp('发电机G1的短路电流IG2=');
disp(IG2);
%线路电流...................................................................
Iij=zeros(n);
for i=1:n
for j=1:n;
Iij(i,j)=(dU(i)-dU(j))*Y(i,j);
end
end
disp('各支路电流矩阵Iij=');
disp(Iij);
%不计负荷，短路前所有电动势，电压均为1，电流为零
% n 为节点数
If3=Uf0(3)/Zf(f);
disp('短路点电流If3=');
disp(9;短路时电压故障分量dU=');
disp(dU);
%运用叠加原理求个节点残余电压...................................................
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)-Y(p,q);
Y(q,q)=Y(q,q)-Y(p,q);
end
%节点导纳矩阵...............................................................
disp('短路分量网络的导纳矩阵Y=');
U=Uf0+dU;
disp('各节点残余电压U=');
disp(U);
%发电机电流................................................................
IG1=(-dU(1))/X(1,2);
disp('发电机G1的短路电流IG1=');
disp(Y);
%节点阻抗矩阵...............................................................
Z=Y^(-1);
disp('短路分量网络的阻抗矩阵Z=');
disp(Z);
%故障点f的选定，及其对应故障分量网络的节点注入电流矩阵I.......................
1,2,0.163+0.754i,0];
% X 为发电机节点号阻抗参数矩阵
% X(i,1)为节点号;
% X(i,2)为发电机电抗;
X=[1,0.15i;
2,0.075i;
3,0;
4,0;
5,0];
%短路前各节点电压Uf0,不考虑负荷
Uf0=[1;1;1;1;1];