1.算法仿真效果
matlab2022a仿真结果如下：
2.算法涉及理论知识概要
在电力工程中，“潮流”还特指电网各处电压（包括幅值与相角）、有功功率、无功功率等的分布。潮流的分布是运行调度单位和维修部门所必须知道的事项。
而潮流计算，是指给定电网中一些参数、已知值和未知值中假设的初始值，通过重复迭代，最终求出潮流分布的精确值，常用方法有牛顿-拉夫逊法和PQ分解法。
在发电机母线上功率被注入网络；而在变（配）电站上接入负荷；其间，功率在网络中流动。对于这种流动的功率，电力生产部门称为潮流（POWER FLOW)。
潮流：电力系统中 电压（各节点）、功率（有功、无功）（各支路）的稳态分布潮流计算---电力系统分析中的一种最基本的计算，根据给定的运行参数确定系统的运行状态，如计算网络中各节点的电压（幅值和相角）和各支路中的功率分布及损耗。
电力系统潮流计算是电力系统最基本的计算，也是最重要的计算。所谓潮流计算，就是已知电网的接线方式与参数及运行条件，计算电力系统稳态运行各母线电压、各支路电流、功率及网损。对于正在运行的电力系统，通过潮流计算可以判断电网母线电压、支路电流和功率是否越限，如果有越限，就应采取措施，调整运行方式。对于正在规划的电力系统，通过潮流计算，可以为选择电网供电方案和电气设备提供依据。潮流计算还可以为继电保护和自动装置整定计算、电力系统故障计算和稳定计算等提供原始数据。
表征电力系统运行状态的参量。包括电力系统中各节点和支路中的电压、电流和功率的流向及分布。在实用上，一般是指稳态运行方式下的静态潮流。合理的潮流分布是电力系统运行的基本要求,其要点为:①运行中的各种电工设备所承受的电压应保持在允许范围内，各种元件所通过的电流应不超过其额定电流，以保证设备和元件的安全；②应尽量使全网的损耗最小，达到经济运行的目的；③正常运行的电力系统应满足静态稳定和暂态稳定的要求。并有一定的稳定储备，不发生异常振荡现象。为此就要求电力系统运行调度人员随时密切监视并调整潮流分布。现代电力系统潮流分布的监视和调整是通过以在线计算机为中心的调度自动化系统来实现的。电力系统潮流的计算和分析是电力系统运行和规划工作的基础。运行中的电力系统，通过潮流计算可以预知,随着各种电源和负荷的变化以及网络结构的改变,网络所有母线的电压是否能保持在允许范围内，各种元件是否会出现过负荷而危及系统的安全，从而进一步研究和制订相应的安全措施。规划中的电力系统，通过潮流计算，可以检验所提出的网络规划方案能否满足各种运行方式的要求，以便制定出既满足未来供电负荷增长的需求，又保证安全稳定运行的网络规划方案。
潮流计算的一般提法是：已知电力网络的结构和参数，已知各负荷点、电源点吸取或发出的有功功率和无功功率（PQ节点），给定电压控制点的电压幅值和有功功率（PV节点），对指定的一个平衡节点给定其电压幅值和相位角（Vθ点）,求解全网各节点电压幅值和相位角，并进一步算出各支路的功率分布和网络损耗。求解潮流问题的基本方程式是节点功率平衡方程。若全网有n个节点，对其中任一节点,可写出其节点功率平衡方程式i=1,2,…,n式中Pi、Qi分别为节点注入有功功率和无功功率,妭i为节点电压相量，Yik为节点导纳矩阵元素。这一方程描述了节点电压同功率之间的非线性关系，是潮流计算的基本方程式。对潮流计算的数字计算机求解方法提出的基本要求是：①计算速度快；②占用存储量少；③收敛性好；④方法简单。
数值解法 潮流计算在数学上是求解一组非线性方程，基本的方法是迭代法。首先发展的潮流问题数字解法是导纳矩阵迭代法。它占用计算机存储量少，适合于计算机发展初期阶段的实际条件，其缺点是收敛性较差。其后发展了阻抗矩阵迭代法，克服了导纳矩阵迭代法收敛性差的缺点，但对大电力系统的计算，占用计算机存储量大。
3.MATLAB核心程序
%求差.
b=xlsread('TEST2.xlsx','C2:F31');
precision=1; %误差
t=0; %迭代次数
pq=24; %pq节点数
%pq=4;
U=b(pq+2:N1,1);
%开始牛拉法
while precision>0.00001
P1=zeros(pq,1);
Q=zeros(pq,1);
P2=zeros(N1-pq-1,1);
U2=zeros(N1-pq-1,1);
deltpqu=zeros(2*(N1-1),1);
deltu2=zeros(N1-pq-1,1);
for i=1:pq
for j=1:N1
.....................................................
end
deltp1=b(2:pq+1,3)-P1(1:pq,1);
deltq =b(2:pq+1,4)-Q(1:pq,1);
deltp2=b(pq+2:N1,3)-P2(1:N1-pq-1,1);
for i=1:N1-pq-1
deltu2(i,1)=U(i,1)^2-U2(i,1);
end
deltu2;
for i=1:pq
deltpqu(2*i-1,1)=deltp1(i);
deltpqu(2*i,1)=deltq(i);
end
for i=1:N1-pq-1
deltpqu(2*pq+2*i-1,1)=deltp2(i);
deltpqu(2*pq+2*i,1)=deltu2(i);
end
deltpqu; %delta
%求雅克比
aii=zeros(N1);bii=zeros(N1);
H1=zeros(N1);N11=zeros(N1);
for i=1:N1
for j=1:N1
aii(i,1)=aii(i,1)+G(i,j)*b(j,1)-B(i,j)*b(j,2);
bii(i,1)=bii(i,1)+G(i,j)*b(j,2)+B(i,j)*b(j,1);
end
end
%H、N
...................................................
Jacobi;
%修正
Correction=Jacobi^-1*deltpqu;
for i=1:N1-1
b(i+1,1)=b(i+1,1)+Correction(2*i);
b(i+1,2)=b(i+1,2)+Correction(2*i-1);
end
precision=max(abs(deltpqu));
%precision=max(abs(Correction));
t=t+1;
b;
end
%平衡节点功率
pp=0;
for j=1:N1
pp=pp+conj(Y(1,j))*(b(j,1)-b(j,2)*1i);
end
b(1,3)=real(b(1,1)*pp);
b(1,4)=imag(b(1,1)*pp);
%各线路功率
S=zeros(N1); y=zeros(N1); y0=zeros(N1);
for i=1:N1
for j=1:N1
if i~=j
y(i,j)=-Y(i,j);
end
end
end
y;
for m1=1:L1
I=B1(m1,1);J=B1(m1,2);R=B1(m1,3);X=B1(m1,4);k=B1(m1,5);
if I*J>0
y0(I,J)=k*1i; y0(J,I)=k*1i;
end
if I*J<0
if I<0
t=I;
I=J;
J=t;
end
J=abs(J);
if k<0
k=-1/k;
end
y0(I,J)=(R/(R^2+X^2)-X/(R^2+X^2)*1i)*(k-1)/k;
y0(J,I)=(R/(R^2+X^2)-X/(R^2+X^2)*1i)*(1-k)/k^2;
end
end