电力系统潮流计算代码

合集下载

电力系统潮流计算matlab程序

电力系统潮流计算matlab程序电力系统潮流计算是电力系统运行和规划中的重要环节，它用于计算电力系统中各节点的电压、功率和电流等参数。

随着电力系统规模的不断扩大和复杂性的增加，传统的手工计算方法已经无法满足需求，因此，利用计算机编程进行潮流计算成为了一种必要的选择。

Matlab是一种功能强大的科学计算软件，它提供了丰富的数学函数和工具箱，可以方便地进行电力系统潮流计算。

下面我将介绍一下如何使用Matlab编写电力系统潮流计算程序。

首先，我们需要建立电力系统的节点模型。

节点模型是电力系统中各节点的电压、功率和电流等参数的数学表示。

在Matlab中，我们可以使用矩阵来表示节点模型。

假设电力系统有n个节点，我们可以定义一个n×n的复数矩阵Y来表示节点之间的导纳关系，其中Y(i,j)表示节点i和节点j之间的导纳。

同时，我们还需要定义一个n×1的复数向量V来表示各节点的电压，其中V(i)表示节点i的电压。

接下来，我们需要编写潮流计算的主程序。

主程序的主要功能是根据节点模型和潮流计算算法，计算出各节点的电压、功率和电流等参数。

在Matlab中，我们可以使用循环语句和矩阵运算来实现潮流计算。

具体的计算过程可以参考电力系统潮流计算的算法。

在编写主程序之前，我们还需要定义一些输入参数，如电力系统的节点数、发电机节点和负荷节点等。

这些参数可以通过用户输入或者读取文件的方式获取。

同时，我们还需要定义一些输出参数，如各节点的电压、功率和电流等。

这些参数可以通过矩阵运算和循环语句计算得到，并输出到文件或者显示在屏幕上。

最后，我们需要进行程序的测试和调试。

可以通过输入一些测试数据，运行程序并检查输出结果是否正确。

如果发现程序有错误或者结果不准确，可以通过调试工具和打印调试信息的方式进行调试。

总之，利用Matlab编写电力系统潮流计算程序可以提高计算效率和准确性，为电力系统的运行和规划提供有力的支持。

当然，编写一个完整的潮流计算程序需要考虑很多细节和特殊情况，这需要有一定的电力系统和编程知识。

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

程序的稳定性分析
程序在不同计算机上的运行结果是否一致。
程序运行过程中，输入数据的变化对输出结果的影响程度。
程序在长时间运行过程中，输出结果是否保持稳定。
程序在处理异常情况时的表现和稳定性。
程序的扩展性分析
代码可读性：C语言程序应具备良好的可读性，方便后续维护和修改算法效率：C语言程序应采用高效的算法，提高计算速度内存占用：C语言程序应合理利用内存，避免内存泄漏和不必要的内存占用扩展性：C语言程序应具备良好的扩展性，方便添加新功能和优化性能
THANK YOU
汇报人：XX
程序的异常处理说明
异常类型：输入错误、计算错误、内存不足等
异常处理方式：使用try-catch 语句进行异常捕获和处理
异常处理流程：当异常发生时，程序会输出错误信息并终止运行
异常处理结果：确保程序在遇到异常时能够正确处理并给出相应的提示信息
C语言程序应用示例
示例程序的输入数据格式
添加标题
添加标题
添加标题Βιβλιοθήκη 输入输出函数：用于数据的输入和输出
函数：可重复使用的代码块，具有特定的功能
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;

电力系统潮流计算完整程序及详细理论说明

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

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

此程序是一个通用的程序，只需要修改输入数据的子函数（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。

电力系统分析潮流计算代码

h0 (n) = u (m)*u (n)*(G(m,n)*sin(delt (m)-delt (n))-B(m,n)*cos(delt (m)-delt (n)));
n0 (n) =-u (m)*u (n)*(G(m,n)*cos(delt (m)-delt (n))+B(m,n)*sin(delt (m)-delt (n)));
end
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
precision=max3;
end
end
end
k-1,delt',u'
%the following program is used to calculate the S5 and Smn
for n=1:N1+1
U (n) =u (n)* (cos (delt(n))+j*sin (delt (n)));
for n=1:N1+1
if m==n
H(m,m)=u(m)^2*B(m,m)+qt(m);
N(m,m)=-u(m)^2*G(m,m)-pt(m);
J(m,m)=u(m)^2*G(m,m)-pt(m);
L(m,m)=u(m)^2*B(m,m)-qt(m);
G(2,4)=0; B(2,4)=0;
G(2,5)=-1.250;B(2,5)=3.750;
G(3,1)=-1.667;B(3,1)=5.000;

电力系统分析潮流计算matlab

目录：一、软件需求说明书......................................................... .. (3)二、概要设计说明书......................................................... .. (4)1、编写潮流计算程序......................................................... . (4)2、数据的输入测试......................................................... .. (4)3、运行得出结果......................................................... (4)4、进行实验结果验证......................................................... . (4)三、详细设计说明书......................................................... .. (5)1、数据导入模块......................................................... (5)2、节点导纳矩阵模块......................................................... . (5)3、编号判断模块......................................................... (5)4、收敛条件判定模块......................................................... .. (5)5、雅可比矩阵模块......................................................... (5)6、迭代计算模块......................................................... . (5)7、计算输出参数模块......................................................... .. (5)四、程序代码......................................................... .. (6)五、最测试例......................................................... (15)1、输入结果......................................................... (15)2、输出结果......................................................... (15)3、结果验证......................................................... (15)一、软件需求说明书本次设计利用MATLAB/C++/C（使用MATLAB）编程工具编写潮流计算，实现对节点电压和功率分布的求取。

Matlab实现潮流计算程序

程序代码如下:111111.%读入数据clcclearfilename='123.txt';a=textread(filename)n=a(1,1);pinghengjd=a(1,2);phjddianya=a(1,3);jingdu=a(1,4);b=zeros(1,9);j1=0;[m1,n1]=size(a);for i1=1:m1if a(i1,1)==0j1=j1+1;b(j1)=i1;endendb;%矩阵分块a1=a(b(1)+1:b(2)-b(1)+1,1:n1);a2=a(b(2)+1:b(3)-1,1:n1);a3=a(b(3)+1:b(4)-1,1:n1);a4=a(b(4)+1:b(5)-1,1:n1);a5=a(b(5)+1:b(6)-1,1:n1);%设置初值vcz=1;dcz=0;kmax=20;k1=0;%求节点导纳矩阵a11=zeros(4,6);for i0=1:3for j0=1:6a11(i0,j0)=a1(i0,j0);a11(4,j0)=a2(1,j0);endenda11;linei=a11(1:4,2);linej=a11(1:4,3);liner=a11(1:4,4);linex=a11(1:4,5);lineb=a11(1:4,6);branchi=0;branchj=0;branchb=0;G=zeros(4,4);B=zeros(4,4);for k=1:4i2=linei(k,1);j2=linej(k,1);r=liner(k,1);x=linex(k,1);b=0;GIJ=r/(r*r+x*x);BIJ=-x/(r*r+x*x);if k>=4 & lineb(k)~=0k0=lineb(k);G(i2,j2)=-GIJ/k0;G(j2,i2)=G(i2,j2);B(i2,j2)=-BIJ/k0;B(j2,i2)=B(i2,j2);G(i2,i2)=G(i2,i2)+GIJ/k0/k0; B(i2,i2)=B(i2,i2)+BIJ/k0/k0;elseG(j2,i2)=-GIJ;G(i2,j2)=G(j2,i2);B(j2,i2)=-BIJ;B(i2,j2)=B(j2,i2);G(i2,i2)=G(i2,i2)+GIJ;b=lineb(k);B(i2,i2)=B(i2,i2)+BIJ+b;endG(j2,j2)=G(j2,j2)+GIJ;B(j2,j2)=B(j2,j2)+BIJ+b;endG;B;B=B.*i;Yf=G+BY=abs(Yf);alf=angle(Yf);%赋Jacobian矩阵参数P=zeros(n,1);Q=zeros(n,1);Pd=zeros(1,n);Qd=zeros(1,n);dP=zeros(1,n);dQ=zeros(1,n);PG=a4(:,3);PD=a4(:,5);QG=a4(:,4);QD=a4(:,6);i8=a4(:,2);for j8=1:length(i8)P(i8(j8))=PG(i8(j8))-PD(i8(j8));Q(i8(j8))=QG(i8(j8))-QD(i8(j8));enddelt=zeros(n,1);V=ones(n,1);V(3)=1.10;V(4)=1.05;ddelt=zeros(n,1);dV=zeros(n,1);A=zeros(2*n,2*n);B=zeros(2*n,1);Jacobian=Jaco(V,delt,n,Y,alf)%求取矩阵功率for j5=1:kmaxdisp(['第' int2str(j5) '次计算结果'])if k>=kmaxbreakendfor i10=1:4Pd(i10)=0;Qd(i10)=0;for j10=1:nPd(i10)=Pd(i10)+V(i10)*Y(i10,j10)*V(j10)*cos(d elt(i10)-delt(j10)-alf(i10,j10));Qd(i10)=Qd(i10)+V(i10)*Y(i10,j10)*V(j10)*sin(d elt(i10)-delt(j10)-alf(i10,j10));endendfor i4=1:3dP(i4)=P(i4)-Pd(i4);endfor j4=1:2dQ(j4)=Q(j4)-Qd(j4);endA=Jaco(V,delt,n,Y,alf)for i14=1:nB(i14*2-1)=-dP(i14);B(i14*2)=-dQ(i14);endif max(abs(B))>jingduX=A\B;for i16=1:nddelt(i16)=X(2*i16-1);dV(i16)=X(2*i16)*V(i16);endV=V+dVdelt=delt+ddeltelsebreakenddisp('----------------')end%流氓算法% for ii=1:2% V(ii)=V(ii)+dV(ii);% end% V222222.function A=Jaco(V,delt,n,Y,alf)%计算Jacobian矩阵for i7=1:nHd1(i7)=0;Jd1(i7)=0;for j7=1:nHd1(i7)=Hd1(i7)+V(i7)*Y(i7,j7)*V(j7)*sin(delt(i7)-delt(j7)-alf(i7,j7));Jd1(i7)=Jd1(i7)+V(i7)*Y(i7,j7)*V(j7)*cos(delt(i7)-delt(j7)-alf(i7,j7));endendfor i6=1:nfor j6=1:nif i6~=j6H(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*sin(delt(i6)-delt(j6)-alf(i6,j6));N(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6));J(i6,j6)=-N(i6,j6);L(i6,j6)=H(i6,j6);elseH(i6,i6)=Hd1(i6)-V(i6)*Y(i6,i6)*V(i6)*sin(delt(i6)-delt(j6)-alf(i6,j6));J(i6,j6)=-Jd1(i6)+V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6));N(i6,j6)=-Jd1(i6)-V(i6)*Y(i6,i6)*V(i6)*cos(alf(i6,i6));L(i6,i6)=-Hd1(i6)+V(i6)*Y(i6,i6)*V(i6)*sin(alf(i6,i6));endendend%修正Jacobian矩阵for j9=3for i9=1:nN(i9,j9)=0;L(i9,j9)=0;J(j9,i9)=0;L(j9,i9)=0;endendL(j9,j9)=1;for j9=4for i9=1:nH(i9,j9)=0;N(i9,j9)=0;J(i9,j9)=0;L(i9,j9)=0;H(j9,i9)=0;N(j9,i9)=0;J(j9,i9)=0;L(j9,i9)=0;endendH(j9,j9)=1;L(j9,j9)=1;%Jaco=[H N;J L];%Jaco=zeros(2*n,2*n);for i11=1:nfor j11=1:nJaco(2*i11-1,2*j11-1)=H(i11,j11); Jaco(2*i11-1,2*j11)=N(i11,j11); Jaco(2*i11,2*j11-1)=J(i11,j11);Jaco(2*i11,2*j11)=L(i11,j11);endendA=Jaco;33333.数据:4 4 1.05 0.000011 12 0.1 0.40 0.015282 1 4 0.12 0.50 0.019203 24 0.08 0.40 0.014131 1 3 0 0.3 0.909090911 1 0 0 0.30 0.182 2 0 0 0.55 0.133 3 0.5 0 0 01 3 1.10 0 0。

电力系统分析潮流计算C语言编程-pq分解法2

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");

电力系统潮流计算的MATLAB辅助程序设计-潮流计算程序

电力系统潮流计算的MATLAB辅助程序设计潮流计算，通常指负荷潮流，是电力系统分析和设计的主要组成部分，对系统规划、安全运行、经济调度和电力公司的功率交换非常重要。

此外，潮流计算还是其它电力系统分析的基础，比如暂态稳定，突发事件处理等。

现代电力系统潮流计算的方法主要:高斯法、牛顿法、快速解耦法和MATLAB的M语言编写的MATPOWER4.1，这里主要介绍高斯法、牛顿法和快速解耦法.高斯法的程序是lfgauss，其与lfybus、busout和lineflow程序联合使用求解潮流功率。

lfybus、busout和lineflow程序也可与牛顿法的lfnewton程序和快速解耦法的decouple程序联合使用。

(读者可以到MATPOWER主页下载MATPOWER4.1，然后将其解压到MATLAB目录下，即可使用该软件进行潮流计算）一、高斯—赛德尔法潮流计算使用的程序：高斯—赛德法的具体使用方法读者可参考后面的实例,这里仅介绍各程序的编写格式：lfgauss：该程序是用高斯法对实际电力系统进行潮流计算，需要用到busdata和linedata两个文件。

程序设计为输入负荷和发电机的有功MW和无功Mvar，以及节点电压标幺值和相角的角度值。

根据所选复功率为基准值将负荷和发电机的功率转换为标幺值。

对于PV节点，如发电机节点,要提供一个无功功率限定值。

当给定电压过高或过低时，无功功率可能超出功率限定值。

在几次迭代之后(高斯—塞德尔迭代为10次），需要检查一次发电机节点的无功出力，如果接近限定值,电压幅值进行上下5%的调整,使得无功保持在限定值内。

lfybus：这个程序需要输入线路参数、变压器参数以及变压器分接头参数。

并将这些参数放在名为linedata的文件中。

这个程序将阻抗转换为导纳,并得到节点导纳矩阵.busout：该程序以表格形式输出结果,节点输出包括电压幅值和相角，发电机和负荷的有功和无功功率，以及并联电容器或电抗器的有功和无功功率。

电力系统潮流分析matlab代码

functiontisco% 这是一个电力系统潮流计算的程序n=input( '\n 请输入节点数:n=') ;m =input( '\n 请输入支路数: m =') ;ph=input( '\n 请输入平衡母线的节点号:ph=') ;B1=input( '\n 请输入支路信息: B1=') ;% 它以矩阵形式存贮支路的情况, 每行存贮一条支路% 第一列存贮支路的一个端点% 第二列存贮支路的另一个端点% 第三列存贮支路的阻抗% 第四列存贮支路的对地导纳% 第五列存贮变压器的变比, 注意支路为1% 第六列存贮支路的序号B2=input( '\n 请输入节点信息: B2=') ;% 第一列为电源侧的功率% 第二列为负荷侧的功率% 第三列为该点的电压值% 第四列为该点的类型: 1 为PQ 节点,2 为PV 节点,3 为平衡节点A =input( '\n 请输入节点号及对地阻抗: A =') ;ip=input( '\n 请输入修正值: ip=') ;% ip为修正值Y =zeros( n) ;e=zeros( 1,n) ;f=zeros( 1,n) ;no=2*ph-1;fori=1:nif A( i,2) ~=0p=A( i,1) ;Y(p,p) =1./A( i,2) ;endendfori=1:mp=B1( i,1) ;q=B1( i,2) ;Y( p,p) =Y( p,p) +1./( B1( i,3) *B1( i,5) ^2+B1( i,4) ./2;Y( 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,4) ./2;endG =real( Y ) ;B=imag( Y ) ;fori=1:ne( i) =real( B2( i,3) ) ;f( i) =imag( B2( i,3) ) ;S(i) =B2( i,1) -B2( i,2) ;V(i) =B2( i,3) ;endP=real( S) ;Q =imag( S) ;[ C ,D ,D F ] =xxf( G ,B ,e,f,P ,Q ,n,B2,ph,V ,no) ;J=jacci( Y ,G ,B ,P ,Q ,e,f,V ,C ,D ,B2,n,ph,no) ; [ De,D f] =hxf( J,D F,ph,n,no) ;t=0;whilem ax( abs( D e) ) >ip& m ax( abs( D f) ) >ipt=t+1;e=e+D e;f=f+D f;[ C ,D ,D F] =xxf( G ,B ,e,f,P,Q ,n,B2,ph,V ,no) ;J=jacci( Y ,G ,B ,P,Q ,e,f,V ,C ,D ,B2,n,ph,no) ; [ De,D f] =hxf( J,D F ,ph,n,no) ;endv=e'+f'*j;fori=1:nhh( i) =conj( Y( ph,i) *v( i) ) ;endS(ph) =sum( hh) *v( ph) ;B2( ph,1) =S( ph) ;V =abs( v) ;jd=angle( v) *180/pi;result1=[ A( :,1) ,real( v) ,imag( v) ,V ,jd,real( S.') , imag( S.') ,real( B2( :,1) ) ,imag( B2( :,1) ) ,real ( B2( :,2) ) ,imag( B2( :,2) ) ] ;fori=1:ma( i) =conj( ( v( B1( i,1) ) /B1( i,5) -v( B1( i,2) ) ) /B1( i,3) ) ;b( i) =v( B1( i,1) ) *a( i) -j*B1( i,4) *v( B1( i,1) ) ^2/2;c( i) =-v( B1( i,2) ) *a( i) -j*B1( i,4) *v( B1( i,2) ) ^2/2;endresult2=[ B1( :,6) ,B1( :,1) ,B1( :,2) ,real( b.') ,imag ( b.') ,real( c.') ,imag( c.') ,real( b.'+c.') ,imag( b.'+c.') ] ;printout( t,result1,S,b,c,result2) ;typeresult.mfunction[ C,D ,D F] =xxf( G ,B,e,f,P,Q ,n,B2,ph,V ,no) % 该子程序是用来求取D Ffori=1:nifi~=phC(i) =0;D(i) =0;for j=1:nC( i) =C( i) +G( i,j) *e( j) -B( i,j) *f( j) ;D( i) =D( i) +G( i,j) *f( j) +B( i,j) *e( j) ;endP1=C(i) *e( i) +D( i) *f( i) ;Q 1=C(i) *f( i) -D( i) *e( i) ;V 1=e(i) ^2+f( i) ^2;if B2( i,4) ==2p=2*i-1;D F( p) =P( i) -P1;p=p+1;D F( p) =V( i) ^2-V 1^2;elsep=2*i-1;D F( p) =P( i) -P1;p=p+1;D F( p) =Q( i) -Q 1;endendendD F=D F';ifph~=nD F( no,:) =[ ] ;D F( no,:) =[ ] ;endfunction [ D e,D f] =hxf( J,D F,ph,n,no)% 该子函数是为求取D e D fD X =J\D F;D X 1=D X ;x1=length( D X 1) ;ifph~=nD X( no) =0;D X( no+1) =0;fori=( no+2) :( x1+2)D X(i) =D X 1( i-2) ;endelseD X =[ D X 1,0,0] ;endk=0;[ x,y] =size( D X ) ;fori=1:2:xk=k+1;D f( k) =D X( i) ;D e( k) =D X( i+1) ;endfunctionJ=jacci( Y ,G ,B,P,Q ,e,f,V ,C,D ,B2,n,ph,no) % 该子程序是用来求取jacci矩阵fori=1:nswitch B2( i,4)case 3continuecase 1for j=1:nif j- =i& j- =phX 1=G( i,j) *f( i) -B( i,j) *e( i) ;X 2=G( i,j) *e( i) +B( i,j) *f( i) ;X 3=-X 2;X 4=X 1;p=2*i-1;q=2*j-1;J(p,q) =X 1;m =p+1;J( m ,q) =X 3;q=q+1;J(p,q) =X 2;J( m ,q) =X 4;else if j==i& j~=phX 1=D( i) +G( i,i) *f( i) -B( i,i) *e( i) ;X 2=C( i) +G( i,i) *e( i) +B( i,i) *f( i) ;X 3=C( i) -G( i,i) *e( i) -B( i,i) *f( i) ;X 4=-D( i) +G( i,i) *f( i) -B( i,i) *e( i) ;p=2*i-1;q=2*j-1;J(p,q) =X 1;m =p+1;J( m ,q) =X 3;q=q+1;J(p,q) =X 2;J( m ,q) =X 4;endendendcase 2for j=1:nif j~=i& j~=phX 1=G( i,j) *f( i) -B( i,j) *e( i) ;X 2=G( i,j) *e( i) +B( i,j) *f( i) ;X 3=0;X 4=0;p=2*i-1;q=2*j-1;J(p,q) =X 1;m =p+1;J( m ,q) =X 3;q=q+1;J(p,q) =X 2;J( m ,q) =X 4;else if j==i& j~=phX 1=D( i) +G( i,i) *f( i) -B( i,i) *e( i) ; X 2=C( i) +G( i,i) *e( i) +B( i,i) *f( i) ; X 3=0;X 4=0;p=2*i-1;q=2*j-1;J(p,q) =X 1;m =p+1;J( m ,q) =X 3;q=q+1;J(p,q) =X 2;J( m ,q) =X 4;endendendendendifph~=nJ( no,:) =[ ] ;J( no,:) =[ ] ;J( :,no) =[ ] ;J( :,no) =[ ] ;end。

电力系统分析潮流计算最终完整版

电力系统分析潮流计算实验报告姓名：XXXXXX 学号：XXXXXXXXXX 班级：XXXXXXXX一、实验目的掌握潮流计算计算机算法的方法，熟悉MATLAB的程序调试方法。

二、实验准备根据课程内容，熟悉MATLAB软件的使用方法，自行学习MATLAB程序的基础语法，并根据所学知识编写潮流计算牛顿拉夫逊法(或PQ分解法) 的计算程序，用相应的算例在MATLAB上进行计算、调试和验证。

三、实验要求每人一组，在实验课时内，调试和修改运行程序，用算例计算输出潮流结果。

四、程序流程五、实验程序%本程序的功能是用牛拉法进行潮流计算%原理介绍详见鞠平著《电气工程》%默认数据为鞠平著《电气工程》例8.4所示数据%B1是支路参数矩阵%第一列和第二列是节点编号。

节点编号由小到大编写%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号%第三列为支路的串列阻抗参数，含变压器支路此值为变压器短路电抗%第四列为支路的对地导纳参数，含变压器支路此值不代入计算%第五烈为含变压器支路的变压器的变比，变压器非标准电压比%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，“0”为不含有变压器%B2为节点参数矩阵%第一列为节点注入发电功率参数%第二列为节点负荷功率参数%第三列为节点电压参数%第四列%第五列%第六列为节点类型参数，“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数%X为节点号和对地参数矩阵%第一列为节点编号%第二列为节点对地参数%默认算例% n=4;% n1=4;% isb=4;% pr=0.00001;% B1=[1 2 0.1667i 0 0.8864 1;1 3 0.1302+0.2479i 0.0258i 1 0;1 4 0.1736+0.3306i 0.0344i 1 0;3 4 0.2603+0.4959i 0.0518i 1 0];% B2=[0 0 1 0 0 2;0 -0.5-0.3i 1 0 0 2;0.2 0 1.05 0 0 3;0 -0.15-0.1i 1.05 0 0 1];% X=[1 0;2 0.05i;3 0;4 0];clear;clc;num=input('是否采用默认数据？(1-默认数据;2-手动输入)');if num==1n=4;n1=4;isb=4;pr=0.00001;B1=[1 2 0.1667i 0 0.8864 1;1 3 0.1302+0.2479i 0.0258i 1 0;1 4 0.1736+0.3306i 0.0344i 1 0;3 4 0.2603+0.4959i 0.0518i 1 0];B2=[0 0 1 0 0 2;0 -0.5-0.3i 1 0 0 2;0.2 0 1.05 0 0 3;0 -0.15-0.1i 1.05 0 0 1];X=[1 0;2 0.05i;3 0;4 0];elsen=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');X=input('节点号和对地参数:X=');endTimes=1; %迭代次数%创建节点导纳矩阵Y=zeros(n);for i=1:n1if B1(i,6)==0 %不含变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else %含有变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-B1(i,5)/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+B1(i,5)/B1(i,3)+(1-B1(i,5))/B1(i,3);Y(q,q)=Y(q,q)+B1(i,5)/B1(i,3)+(B1(i,5)*(B1(i,5)-1))/B1(i,3);endendfor i=1:n1Y(i,i)=Y(i,i)+X(i,2); %计及补偿电容电纳enddisp('导纳矩阵为：');disp(Y); %显示导纳矩阵%初始化OrgS、DetaSOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);%创建OrgS，用于存储初始功率参数h=0;j=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2 %不是平衡点&是PQ点h=h+1;for j=1:n%公式8-74%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej)%Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real (Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j ,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:n %对PV节点的处理，注意这时不可再将h初始化为0if i~=isb&B2(i,6)==3 %不是平衡点&是PV点h=h+1;for j=1:n%公式8-75-a%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej)%Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real (Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j ))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendend%创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:nif B2(i,6)==3t=t+1;PVU(t,1)=B2(i,3);endend%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1); %delPiDetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1); %delQiendendt=0;for i=1:n %对PV节点的处理，注意这时不可再将h初始化为0if i~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1); %delPiDetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2; %delUi endend% DetaS%创建I，用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));%conj求共轭endend%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理if B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+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));else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));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:n %对PV节点的处理if B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendenddisp('初始雅可比矩阵为：');disp(Jacbi);%求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS; %inv矩阵求逆% DetaU%修正节点电压j=0;for i=1:n %对PQ节点处理if B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:n %对PV节点的处理if B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endend% B2%开始循环**********************************************************************while abs(max(DetaU))>prOrgS=zeros(2*n-2,1);h=0;j=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real (Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j ))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:nif i~=isb&B2(i,6)==3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real (Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j ))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendend% OrgS%创建DetaSh=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);endendt=0;for i=1:nif i~=isb&B2(i,6)==3h=h+1;t=t+1;% DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;endend% DetaS%创建Ii=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));endend% I%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+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)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));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,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendend% JacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;% DetaU%修正节点电压j=0;for i=1:nif B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:nif B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endend% B2Times=Times+1; %迭代次数加1enddisp('迭代次数为：');disp(Times);disp('收敛时电压修正量为：：');disp(DetaU);for k=1:nE(k)=B2(k,3);e(k)=real(E(k));f(k)=imag(E(k));V(k)=sqrt(e(k)^2+f(k)^2);sida(k)=atan(f(k)./e(k))*180./pi;end%=============== 计算各输出量=========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：'); disp(E); %显示各节点的实际电压标幺值E用复数表示disp('-----------------------------------------------------')disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V); %显示各节点的电压大小V的模值disp('-----------------------------------------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida); %显示各节点的电压相for p=1:nfor q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q)); %计算各节点的注入电流的共轭值endS(p)=E(p)*C(p); %计算各节点的功率S = 电压X 注入电流的共轭值enddisp('各节点的功率S为(节点号从小到大排列)：');disp(S); %显示各节点的注入功率Sline=zeros(n1,5);disp('-----------------------------------------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:n1p=B1(i,1);q=B1(i,2);Sline(i,1)=B1(i,1);Sline(i,2)=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))*((1-B1(i,5))/B1(i,3))+(conj(E(p))-conj(E(q)))*(B1(i,5)/B1(i,3)));Siz(i)=Si(p,q);endSSi(p,q)=Si(p,q);Sline(i,3)=Siz(i);ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);enddisp('-----------------------------------------------------');disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:n1p=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))*((B1(i,5)*(B1(i,5)-1))/B1(i,3))+(conj(E(q))-conj(E(p)))*(B1(i,5)/B1(i,3)));Sjy(i)=Sj(q,p);endSSj(q,p)=Sj(q,p);Sline(i,4)=Sjy(i);ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);enddisp('-----------------------------------------------------');disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:n1p=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);DDS(i)=DS(i);Sline(i,5)=DS(i);ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);enddisp('-----------------------------------------------------');disp('各支路首端编号末端编号首端功率末端功率线路损耗');disp(Sline);六、运行结果及其分析是否采用默认数据？(1-默认数据;2-手动输入)1导纳矩阵为：2.9056 -11.5015i 0.0000 + 5.3173i -1.6606 +3.1617i -1.2450 + 2.3710i0.0000 + 5.3173i 0.0000 - 4.6633i 0.0000 + 0.0000i 0.0000 + 0.0000i-1.6606 + 3.1617i 0.0000 + 0.0000i 2.4904 - 4.7039i -0.8298 + 1.5809i-1.2450 + 2.3710i 0.0000 + 0.0000i -0.8298 + 1.5809i 2.0749 - 3.9089i初始雅可比矩阵为：11.1267 2.7603 -5.3173 0 -3.1617 -1.6606-3.0509 11.8762 0 -5.3173 1.6606 -3.1617-5.3173 0 5.3173 0 0 00 -5.3173 0 4.0092 0 0-3.3198 -1.7436 0 0 4.8217 2.69800 0 0 0 0 2.1000迭代次数为：4收敛时电压修正量为：：1.0e-05 *0.0349-0.2445-0.0101-0.5713-0.0931-0.0073各节点的实际电压标幺值E为(节点号从小到大排列)：0.9673 - 0.0655i 1.0252 - 0.1666i 1.0495 - 0.0337i 1.0500 + 0.0000i -----------------------------------------------------各节点的电压大小V为(节点号从小到大排列)：0.9695 1.0387 1.0500 1.0500-----------------------------------------------------各节点的电压相角sida为(节点号从小到大排列)：-3.8734 -9.2315 -1.8419 0各节点的功率S为(节点号从小到大排列)：-0.0000 + 0.0000i -0.5000 - 0.3000i 0.2000 + 0.1969i 0.3277 + 0.0443i -----------------------------------------------------各条支路的首端功率Si为(顺序同您输入B1时一致)：S(1,2)=-0.5-0.30713iS(1,3)=-0.24266-0.197iS(1,4)=-0.25734-0.11013iS(3,4)=-0.055551+0.0017528i-----------------------------------------------------各条支路的末端功率Sj为(顺序同您输入B1时一致)：S(2,1)=0.5+0.24606iS(3,1)=0.25555+0.1952iS(4,1)=0.2712+0.1014iS(4,3)=0.056496-0.057061i-----------------------------------------------------各条支路的功率损耗DS为(顺序同您输入B1时一致)：DS(1,2)=0-0.06107iDS(1,3)=0.012892-0.0018014iDS(1,4)=0.013863-0.0087295iDS(3,4)=0.00094545-0.055308i-----------------------------------------------------各支路首端编号末端编号首端功率末端功率线路损耗1.0000 + 0.0000i2.0000 + 0.0000i -0.5000 - 0.3071i 0.5000 + 0.2461i 0.0000 - 0.0611i 1.0000 + 0.0000i3.0000 + 0.0000i -0.2427 - 0.1970i 0.2556 + 0.1952i 0.0129 - 0.0018i 1.0000 + 0.0000i4.0000 + 0.0000i -0.2573 - 0.1101i 0.2712 + 0.1014i 0.0139 - 0.0087i3.0000 + 0.0000i4.0000 + 0.0000i -0.0556 + 0.0018i 0.0565 - 0.0571i 0.0009 - 0.0553i七、实验体会及感悟通过这次实验，首先让我对matlab软件有了初步的了解，对它强大的矩阵运算能力有了更深的体会，同时掌握了设置断点和断点调试的一般方法，结合课本上的程序流程图和参考资料上的例子单步跟踪调试，再一次的熟悉了牛顿拉夫逊法潮流计算的一般方法和步骤，对计算机计算潮流计算有了更进一步的认识，在学习潮流计算时，虽然依次学习了节点导纳矩阵，功率方程、雅可比矩阵，但不能将它们联系起来，更不知道其中的原委，通过程序的编写，知道了其中的联系，也知道了每个方程、矩阵在计算中的作用。

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

实验目的根据所给的电力系统，编制潮流计算程序，通过计算机进行调试，最后完成一个切实可行的电力系统计算应用程序.通过自己设计电力系统计算程序使同学对电力系统分析有进一步理解，同时加强计算机实际应用能力的训练。

程序计算原理1、概述应用计算机进行电力系统计算，首先要掌握电力系统相应计算的数学模型；其次是运用合理的计算方法；第三则是选择合适的计算机语言编制计算程序。

建立电力系统计算的相关数学模型，就是建立用于描述电力系统相应计算的有关参数间的相互关系的数学方程式。

该数学模型的建立往往要突出问题的主要方面，即考虑影响问题的主要因素，而忽略一些次要因素，使数学模型既能正确地反映实际问题，又使计算不过于复杂。

运用合理的计算方法，就是要求所选用的计算方法能快速准确地得出正确结果，同时还应要求在解算过程中占用内存少，以利提高计算机的解题规模。

选择合适的语言编写程序，就是首先确定用什么计算机语言来编制程序；其次是作出计算的流程图；第三根据流程图用选择的语言编写计算程序.然后上机调试，直到语法上无错误。

本程序采用C语言进行编程。

所编制的程序难免存在逻辑错误,因此先用一个已知结果的系统作为例题进行计算.用程序计算的结果和已知结果相比较,如果结果相差甚远就要逐步分析程序的计算步骤，查出问题的出处；如果结果比较接近,则逐步分析误差来源;直到结果正确为止。

2、电力系统潮流计算的程序算法潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压(幅值及相角）、网络中的功率分布及功率损耗等.目前计算机潮流计算的方法主要有牛顿—拉夫逊算法和PQ分解法。

牛顿-拉夫逊算法是数学上求解非线形方程组的有效方法,具有较好的收敛性,曾经是潮流计算中应用比较普遍的方法。

PQ快速分解法是从牛顿—拉夫逊算法演变而来的，是将纯数学的牛顿—拉夫逊算法与电力系统具体特点相结合并进行简化与改进而得出的。

PQ快速分解法比牛顿—拉夫逊算法大大提高了计算速度和节省了内存，故而本程序以PQ快速分解法进行潮流计算.1）形成节点导纳矩阵（1)自导纳的形成对节点i其自导纳Y ii是节点i以外的所有节点都接地时节点i对地的总导纳。

电力系统潮流计算的C语言实现

//////////////////////////////////////////////////////////////////////// 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[]:分别存放左、右节点号和电抗器节点号。

电力系统潮流计算的MATLAB辅助程序设计,潮流计算程序(精编文档).doc

【最新整理，下载后即可编辑】电力系统潮流计算的MATLAB辅助程序设计潮流计算，通常指负荷潮流，是电力系统分析和设计的主要组成部分，对系统规划、安全运行、经济调度和电力公司的功率交换非常重要。

此外，潮流计算还是其它电力系统分析的基础，比如暂态稳定，突发事件处理等。

现代电力系统潮流计算的方法主要:高斯法、牛顿法、快速解耦法和MATLAB的M语言编写的MATPOWER4.1，这里主要介绍高斯法、牛顿法和快速解耦法。

高斯法的程序是lfgauss，其与lfybus、busout和lineflow程序联合使用求解潮流功率。

lfybus、busout和lineflow程序也可与牛顿法的lfnewton程序和快速解耦法的decouple程序联合使用。

(读者可以到MATPOWER主页下载MATPOWER4.1,然后将其解压到MATLAB目录下，即可使用该软件进行潮流计算)一、高斯-赛德尔法潮流计算使用的程序：高斯-赛德法的具体使用方法读者可参考后面的实例，这里仅介绍各程序的编写格式:lfgauss：该程序是用高斯法对实际电力系统进行潮流计算，需要用到busdata和linedata两个文件。

程序设计为输入负荷和发电机的有功MW和无功Mvar，以及节点电压标幺值和相角的角度值。

根据所选复功率为基准值将负荷和发电机的功率转换为标幺值。

对于PV节点，如发电机节点，要提供一个无功功率限定值。

当给定电压过高或过低时，无功功率可能超出功率限定值。

在几次迭代之后(高斯-塞德尔迭代为10次)，需要检查一次发电机节点的无功出力，如果接近限定值，电压幅值进行上下5%的调整，使得无功保持在限定值内。

lfybus：这个程序需要输入线路参数、变压器参数以及变压器分接头参数。

并将这些参数放在名为linedata的文件中。

这个程序将阻抗转换为导纳，并得到节点导纳矩阵。

busout：该程序以表格形式输出结果，节点输出包括电压幅值和相角，发电机和负荷的有功和无功功率，以及并联电容器或电抗器的有功和无功功率。

电力系统通用潮流计算C语言程序

＃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。

电力系统分析潮流计算

电力系统分析潮流计算电力系统分析是对电力系统运行状态进行研究、分析和评估的一项重要工作。

其中，潮流计算是电力系统分析的一种重要方法，用于计算电力系统中各节点的电压、功率和电流等参数。

本文将详细介绍电力系统潮流计算的原理、方法和应用。

一、电力系统潮流计算的原理电力系统潮流计算是基于潮流方程的求解，潮流方程是描述电力系统各节点电压和相角之间的关系的一组非线性方程。

潮流方程的基本原理是基于电力系统的等效导纳矩阵和节点电压相位差的关系，通过潮流计算可以得到电力系统各节点的电压和功率等参数。

电力系统潮流方程的一般形式如下：\begin{align*}P_i &= \sum_{j=1}^{n}(V_iV_j(G_{ij}\cos(\theta_i-\theta_j)+B_{ij}\sin(\theta_i-\theta_j))) \\Q_i &= \sum_{j=1}^{n}(V_iV_j(G_{ij}\sin(\theta_i-\theta_j)-B_{ij}\cos(\theta_i-\theta_j)))\end{align*}其中，$n$为节点数，$P_i$和$Q_i$表示第i个节点的有功功率和无功功率。

$V_i$和$\theta_i$表示第i个节点的电压和相角。

$G_{ij}$和$B_{ij}$表示节点i和节点j之间的等效导纳。

二、电力系统潮流计算的方法电力系统潮流计算的方法主要包括直接法、迭代法和牛顿-拉夫逊法等。

1.直接法：直接法是一种适用于小规模电力系统的潮流计算方法，它通过直接求解潮流方程来计算电力系统的潮流。

直接法的计算速度快，但对系统规模有一定的限制。

2.迭代法：迭代法是一种常用的潮流计算方法，通常使用高尔顿法或牛顿法。

迭代法通过迭代求解潮流方程来计算电力系统的潮流。

迭代法相对于直接法来说，可以适用于大规模电力系统，但计算时间较长。

3.牛顿-拉夫逊法：牛顿-拉夫逊法是一种高效的潮流计算方法，它通过求解潮流方程的雅可比矩阵来进行迭代计算，可以有效地提高计算速度。

电力系统潮流计算程序(详细)

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

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

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

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

潮流计算程序，相关的原始数据数据数据输入格式如下：%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。

matlab实验电力系统潮流计算

matlab实验电力系统潮流计算电力系统潮流计算是电力系统运行分析的基础，它通过计算电力系统中各节点的电压幅值和相角，以及各支路的功率和电流，来研究电力系统的稳态工作状态。

本文将介绍潮流计算的原理及其在电力系统中的应用。

潮流计算的基本原理是基于电力系统节点间的功率平衡方程，即节点注入功率等于节点负荷消耗功率加上节点发电机注入功率和节点之间传输输电功率的代数和。

潮流计算通常分为直流潮流计算和交流潮流计算两种方法。

直流潮流计算是指忽略电流相位差的计算方法。

在直流潮流计算中，电力系统的节点电压幅值和相角可以用复数来表示，节点注入功率和节点负荷消耗功率也采用复功率的形式。

直流潮流计算的基本方程为：P+iQ = V*(Gcosθ+Bsinθ)其中，P和Q分别表示节点注入有功功率和无功功率，V表示节点电压幅值，θ表示节点电压相角，G和B分别表示节点导纳矩阵的实部和虚部。

交流潮流计算考虑了电流相位差的影响，是更为准确的潮流计算方法。

交流潮流计算通常采用牛顿-拉夫逊法（Newton-Raphson method）进行迭代求解。

该方法以功率不平衡最小为目标，通过迭代计算更新节点电压幅值和相角，直到收敛为止。

潮流计算在电力系统运行和规划中具有重要的应用价值。

首先，潮流计算可以用来评估电力系统的稳态工作状态，包括节点电压幅值和相角、支路功率和电流等信息。

通过分析潮流计算结果，可以发现电力系统中潜在的潮流瓶颈和潮流分布情况，为电网调度和运行提供参考依据。

其次，潮流计算可以用来优化电力系统的运行和规划。

通过分析潮流计算结果，可以确定潮流分布不均衡的节点和支路，进而优化电力系统的输电和变电容量配置，提高电力系统的可靠性、经济性和稳定性。

此外，潮流计算还可用于电力系统的故障分析和重构，对于故障点的电压幅值、相角以及故障后的支路功率和电流进行分析，有助于电力系统故障的定位和恢复。

总的来说，电力系统潮流计算是电力系统运行分析的重要工具，通过计算电力系统中各节点的电压幅值和相角，以及各支路的功率和电流，可以评估电力系统的稳态工作状态，优化电力系统的运行和规划，实现电力系统的安全、稳定和高效运行。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

附录程序的主要代码：n=input('请输入节点数n=');na=input('请输入支路数na=');isb=input('请输入平衡节点母线号isb=');jd=input('请输入误差精度jd=');B1=input('请输入由支路参数形成的矩阵B1=');B2=input('请输入由节点参数形成的矩阵B2=');L=input('请输入由节点号及其对地阻抗形成的矩阵L='); nb=input('请输入P-Q节点数nb=');Y=zeros(n);Z=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n); O=zeros(1,n);for i=1:naif B1(i,6)==0a=B1(i,1);b=B1(i,2);else a=B1(i,2);b=B1(i,1);endY(a,b)=Y(a,b)-1./(B1(i,3)*B1(i,5));Z(a,b)=Z(a,b)-1./(B1(i,3));Y(b,a)=Y(a,b);Z(b,a)=Z(a,b);Y(b,b)=Y(b,b)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;Z(b,b)=Z(b,b)+1./(B1(i,3));Y(a,a)=Y(a,a)+1./(B1(i,3))+B1(i,4)./2;Z(a,a)=Z(a,a)+1./(B1(i,3));endG=real(Y);B=imag(Z);CI=imag(Y);for i=1:nS(i)=B2(i,1)-B2(i,2);CI(i,i)=CI(i,i)+B2(i,5);endP=real(S);Q=imag(S);for i=1:ne(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4);endfor i=1:nif B2(i,6)==2V(i)=sqrt(e(i)^2+f(i)^2);O(i)=atan(f(i)./e(i));endfor i=2:nif i==nB(i,i)=1./B(i,i);elseIT1=i+1;for j1=IT1:nB(i,j1)=B(i,j1)./B(i,i);endB(i,i)=1./B(i,i);for k=i+1:nfor j1=i+1:nB(k,j1)=B(k,j1)-B(k,i)*B(i,j1);endendendenda=0;b=0;for i=1:nif B2(i,6)==2a=a+1;k=0;for j1=1:nif B2(j1,6)==2k=k+1;A(a,k)=CI(i,j1);endendendendfor i=1:nbif i==naA(i,i)=1./A(i,i);elsek=i+1;for j1=k:nbA(i,j1)=A(i,j1)./A(i,i);endA(i,i)=1./A(i,i);for k=i+1:nbfor j1=i+1:nbA(k,j1)=A(k,j1)-A(k,i)*A(i,j1);endendendNT2=1;NT1=0;kp=1;kq=1;K=1;NCT=0;NT3=1;while NT2~=0|NT3~=0NT2=0;NT3=0;for i=1:nif i~=isbC(i)=0;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+CI(i,k)*sin(O(i)-O(k)));endCP1(i)=P(i)-V(i)*C(i);CP(i)=CP1(i)./V(i);NCT=abs(CP1(i));if NCT>=jdNT2=NT2+1;endendendNp(k)=NT2;if NT2~=0for i=2:nCP(i)=B(i,i)*CP(i);if i~=nIT1=i+1;for k=IT1:nCP(k)=CP(k)-B(k,i)*CP(i);endelsefor LZ=3:iL=i+3-LZ;NC4=L-1;for MZ=2:NC4I=NC4+2-MZ;CP(I)=CP(I)-B(I,L)*CP(L);endendendendfor i=2:nO(i)=O(i)-CP(i);endkq=1;L=0;for i=1:nif B2(i,6)==2C(i)=0;L=L+1;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-CI(i,k)*cos(O(i)-O(k)));endDQ1(i)=Q(i)-V(i)*C(i);DQ(L)=DQ1(i)./V(i);NCT=abs(DQ1(i));if NCT>=jdNT3=NT3+1;endendendelsekp=0;if kq~=0;L=0;for i=1:nif B2(i,6)==2C(i)=0;L=L+1;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-CI(i,k)*cos(O(i)-O(k)));endDQ1(i)=Q(i)-V(i)*C(i);DQ(L)=DQ1(i)./V(i);NCT=abs(DQ1(i));endendendendNq(K)=NT3;if NT3~=0L=0;for i=1:nbDQ(i)=A(i,i)*DQ(i);if i==nbfor LZ=2:iL=i+2-LZ;NC4=L-1;for MZ=1:NC4I=NC4+1-MZ;DQ(I)=DQ(I)-A(I,L)*DQ(L);endendelseIT1=i+1;for k=IT1:nbDQ(k)=DQ(k)-A(k,i)*DQ(i);endendendL=0;for i=1:nif B2(i,6)==2L=L+1;V(i)=V(i)-DQ(L);endendkp=1;K=K+1;elsekq=0;if kp~=0K=K+1;endendfor i=1:nDp(K-1,i)=V(i);endenddisp('迭代次数');disp(K);disp('每次没有达到精度要求的有功功率个数为');disp(Np);disp('每次没有达到精度要求的无功功率个数为');disp(Nq);for k=1:nE(k)=V(k)*cos(O(k))+V(k)*sin(O(k))*j;O(k)=O(k)*180./pi;enddisp('各节点的实际电压标幺值E(节点号从小到大排列):'); disp(E);disp('各节点的电压大小V为(节点号从小到大排列):'); disp(V);disp('各节点的电压相角O为(节点号从小到大排列):'); disp(O);for a=1:nC(a)=0;for b=1:nC(a)=C(a)+conj(Y(a,b))*conj(E(b));endS(a)=E(a)*C(a);enddisp('各节点的功率S为(节点号从小到大排列):');disp(S);disp('各条支路的首端功率Si为(顺序同您输入B1时一样):');for i=1:naif B1(i,6)==0a=B1(i,1);b=B1(i,2);elsea=B1(i,2);b=B1(i,1);endSi(a,b)=E(a)*(conj(E(a))*conj(B1(i,4)./2)+(conj(E(a)*B1(i,5))-conj(E(b)))*conj(1./(B1(i,3)*B1(i,5))));disp(Si(a,b));enddisp('各条支路的末端功率Sj为(顺序同您输入B1时一样):');for i=1:naif B1(i,6)==0a=B1(i,1);b=B1(i,2);elsea=B1(i,2);b=B1(i,1);endSj(b,a)=E(b)*(conj(E(b))*conj(B1(i,4)./2)+(conj(E(b)./B1(i,5))-conj(E(a)))*conj(1./(B1(i,3)*B1(i,5))));disp(Sj(b,a));enddisp('各条支路的功率损耗DS为(顺序同您输入B1时一样)::');for i=1:naif B1(i,6)==0a=B1(i,1);b=B1(i,2);elsea=B1(i,2);b=B1(i,1);endDS(i)=Si(a,b)+Sj(b,a);disp(DS(i));endfor i=1:KCs(i)=i;for j=1:nDp(K,j)=Dp(K-1,j);endenddisp('以下是每次迭代后各节点的电压值（如图所示）');plot(Cs,Dp),xlabel('迭代次数'),ylabel('电压'),title('电压迭代次数曲线');。