✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。
🍎 往期回顾关注个人主页：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信。
🔥 内容介绍
一、引言：潮流计算的核心意义与高斯 - 赛德尔算法的价值
1.1 潮流计算：电力系统分析的 “基石”
潮流计算是电力系统规划、运行与控制的核心基础，其本质是求解电力网络的稳态运行状态 —— 通过已知的负荷、电源参数，计算各母线的电压幅值、相角，以及线路功率流向与损耗。无论是电网扩容规划、新能源接入评估，还是调度运行中的安全校验，都离不开潮流计算的支撑。
1.2 两节点系统：潮流计算的 “入门模型”
两节点电力系统是最简单的电力网络拓扑（含 1 个平衡节点、1 个 PQ 节点），虽结构简单，但能直观体现潮流计算的核心逻辑。本文聚焦两节点系统，采用经典的高斯 - 赛德尔（Gauss-Seidel）迭代法，手把手演示 PQ 节点（未知电压幅值和相角）的求解过程，帮助读者快速掌握潮流计算的核心思路。
1.3 高斯 - 赛德尔算法：迭代法的 “经典代表”
高斯 - 赛德尔算法是潮流计算中最基础的迭代算法，核心优势是原理简单、计算量小、易实现—— 基于节点功率平衡方程，通过逐次迭代更新未知节点电压，直至满足收敛条件。虽收敛速度慢于牛顿 - 拉夫逊法，但对于简单系统（如两节点、三节点），是理解潮流计算本质的最佳切入点。
二、核心原理：两节点系统建模与高斯 - 赛德尔迭代公式
2.1 两节点系统拓扑与参数定义
2.1.1 系统结构
两节点系统由母线 1（平衡节点） 、母线 2（PQ 节点） 及连接线路组成，结构如下：
平衡节点（母线 1）：电压幅值 V₁和相角 θ₁已知（通常设为 V₁=1.0∠0° pu，作为电压参考点），承担平衡系统功率损耗的角色；
PQ 节点（母线 2）：注入有功功率 P₂、无功功率 Q₂已知（由负荷或电源决定），需求解电压幅值 V₂和相角 θ₂；
线路参数：线路阻抗 Z=R+jX（或导纳 Y=G+jB=1/Z，B 为线路电纳，忽略时 Y=G+j0）。
2.1.2 节点导纳矩阵（Y 矩阵）
两节点系统的节点导纳矩阵为 2×2 矩阵，核心是描述节点间的电气连接关系：
Image
Image
⛳️ 运行结果
Image
Image
📣 部分代码
%---------Initial settings-----------
clc;
clear;
close all
Convergence_Tolerance=1e-6; %Newton Raphon Convergence Tolerance
disp(['Convergence Tolerance, Ep:',num2str(Convergence_Tolerance)]);
%-----------------Node, tributary data---------------
% Node data
node = ...
[...
1 3 0 0.065 0 0 1.06 0 0
2 1 1.82 0.64 0 0 1 0 0
];
% Matrix formed by node data: column 1 node number, column 2 node type:1 - PQ 2 - PV 3 - balance
% column 3 load active, column 4 load reactive, column 5 generation active, column 6 generation reactive
% column 7 voltage real, column 8 voltage imaginary
% column 9 node-to-ground conductance
% Branch data
branch = ...
[...
%bus_i %bus_j % R % X % Bi/2 % k
1 2 0 0.034 0 1;
];
% bus_i denotes the first node of a branch.
% bus_j denotes the end node of the branch.
% R and X denote the branch resistance and reactance, respectively, and B denotes the dana.
% k denotes the transformer ratio of the first node to the last node.
n=size(node,1); %Number of nodes
l=size(branch,1); %Number of branches
%%--------------------------------------
%Calculating the nodal derivative matrix
%%--------------------------------------
Y=zeros(n); %Calculating the nodal conductance matrix
for j=1:l
s1=branch(j,1);s2=branch(j,2); %s1 is the first segment node of each branch and s2 is the end node
kij=branch(j,6); %Branch transformer ratio
Y(s1,s2)=-1/( (branch(j,3)+1i*branch(j,4))*kij );
Y(s2,s1)=Y(s1,s2);
Y(s1,s1)=Y(s1,s1)+1/( (branch(j,3)+1i*branch(j,4))*kij^2 )+i*branch(j,5);
Y(s2,s2)=Y(s2,s2)+1/(branch(j,3)+1i*branch(j,4))+i*branch(j,5);
end
for a=1:n %--Calculation of conductivity to earth--
Y(a,a)=Y(a,a)+1i*node(a,9);
end
G=real(Y); %Nodal derivative matrix real part matrix
B=imag(Y); %Nodal derivative matrix imaginary part matrix
%%Newton Raphon method of calculating currents
e=node(:,7); %Initial value of the real part of the voltage
f=zeros(n,1); %Initial value of the imaginary part of the voltage
for a=1:n %Search to balance nodes
if node(a,2)==3
slack=a;
end
end
it=0; %Number of iterations
Convergence=[];
while 1
DP=zeros(n,1); DQ=zeros(n,1); %Uneven measurement
Pi=zeros(n,1); Qi=zeros(n,1); %Nodal injected power calculated by the tidal equation
for a=1:n %Calculating Pi, Qi
if node(a,2)==1 || node(a,2)==3 %Determined to be a PQ node
for b=1:n
Pi(a,1)=Pi(a,1)+e(a,1)*( G(a,b)*e(b,1)-B(a,b)*f(b,1) )+ f(a,1)*( G(a,b)*f(b,1)+B(a,b)*e(b,1) );
Qi(a,1)=Qi(a,1)+f(a,1)*( G(a,b)*e(b,1)-B(a,b)*f(b,1) )- e(a,1)*( G(a,b)*f(b,1)+B(a,b)*e(b,1) );
end
DP(a,1)=( node(a,5)-node(a,3) )-Pi(a,1); %Calculating unevenness
DQ(a,1)=( node(a,6)-node(a,4) )-Qi(a,1);
elseif node(a,2)==2 %Determined to be a PV node
for b=1:n
Pi(a,1)=Pi(a,1)+e(a,1)*( G(a,b)*e(b,1)-B(a,b)*f(b,1) )+ f(a,1)*( G(a,b)*f(b,1)+B(a,b)*e(b,1) );
Qi(a,1)=e(a,1)^2+f(a,1)^2;
end
DP(a,1)=( node(a,5)-node(a,3) )-Pi(a,1);
DQ(a,1)=node(a,7)^2-Qi(a,1); %PV node voltage imbalance
end %Balance nodes regardless of
end
DPQ=zeros(2*n,1); %Integration of active, reactive and voltage inequalities into one matrix
for a=1:n
DPQ(2*a-1,1)=DP(a,1);
DPQ(2*a,1)=DQ(a,1);
end
DPQ([2slack-1,2slack],:)=[]; % The slack node is a balancing node and does not participate in iterations, remove
xV=zeros(n,1); %Node voltage amplitude
xR=zeros(n,1); %Node voltage phase angle
for a=1:n
xV(a,1)=sqrt( e(a,1)^2+f(a,1)^2 );
xR(a,1)=atan(f(a,1)/e(a,1))*180/pi;
end
disp(['Number of iteration:',num2str(it),' Max(|DP|,|DQ|):',num2str(max(abs(DPQ))),' Voltage Magnitude at bus 2:',num2str(xV(2,1)),' Voltage Angle at bus 2:',num2str(xR(2,1))]);
%If the inequality is less than the convergence accuracy, then the convergence condition has been reached and the iteration is skipped
Convergence=[Convergence;max(abs(DPQ))];
if max(abs(DPQ))<=Convergence_Tolerance
break;
end
%If accuracy is achieved, the Jacobi matrix and corrections are calculated downwards and the iteration continues
H=zeros(n,n); % Calculation of DP/De
N=zeros(n,n); % Calculation of DP/Df
K=zeros(n,n); % Calculate DQ/De for the PQ node and DU/De for the PV node
L=zeros(n,n); % Calculate DQ/Df for the PQ node and DU/Df for the PV node
for a=1:n
if node(a,2)==1 %Determined to be a PQ node
for b=1:n
if a~=b
H(a,b)=-( G(a,b)*e(a,1)+B(a,b)*f(a,1) );
N(a,b)=B(a,b)*e(a,1)-G(a,b)*f(a,1) ;
K(a,b)=B(a,b)*e(a,1)-G(a,b)*f(a,1) ;
L(a,b)=G(a,b)*e(a,1)+B(a,b)*f(a,1) ;
elseif a==b
for c=1:n
H(a,b)=H(a,b)-( G(a,c)*e(c,1)-B(a,c)*f(c,1) );
N(a,b)=N(a,b)-( G(a,c)*f(c,1)+B(a,c)*e(c,1) );
K(a,b)=K(a,b)+( G(a,c)*f(c,1)+B(a,c)*e(c,1) );
L(a,b)=L(a,b)-( G(a,c)*e(c,1)-B(a,c)*f(c,1) );
end
H(a,b)=H(a,b)-G(a,a)*e(a,1)-B(a,a)*f(a,1);
N(a,b)=N(a,b)+B(a,a)*e(a,1)-G(a,a)*f(a,1);
K(a,b)=K(a,b)+B(a,a)*e(a,1)-G(a,a)*f(a,1);
L(a,b)=L(a,b)+G(a,a)*e(a,1)+B(a,a)*f(a,1);
end
end
elseif node(a,2)==2 %Determined to be a PV node
for b=1:n
if a~=b
H(a,b)=-( G(a,b)*e(a,1)+B(a,b)*f(a,1) );
N(a,b)=B(a,b)*e(a,1)-G(a,b)*f(a,1) ;
K(a,b)=0 ;
L(a,b)=0 ;
elseif a==b
for c=1:n
H(a,b)=H(a,b)-( G(a,c)*e(c,1)-B(a,c)*f(c,1) );
N(a,b)=N(a,b)-( G(a,c)*f(c,1)+B(a,c)*e(c,1) );
end
H(a,b)=H(a,b)-G(a,a)*e(a,1)-B(a,a)*f(a,1);
N(a,b)=N(a,b)+B(a,a)*e(a,1)-G(a,a)*f(a,1);
K(a,b)=-2*e(a,1);
L(a,b)=-2*f(a,1);
end
end
end %Balance nodes regardless of
end
jacobi=zeros(2*n); %Jacobi Matrix
for a=1:n
for b=1:n
jacobi(2*a-1,2*b-1)=H(a,b);
jacobi(2*a-1,2*b)=N(a,b);
jacobi(2*a,2*b-1)=K(a,b);
jacobi(2*a,2*b)=L(a,b);
end
end
jacobi([2*slack-1,2*slack],:)=[]; jacobi(:,[2*slack-1,2*slack])=[]; % The slack node is a balancing node and does not participate in iterations, remove
delt=-jacobi\DPQ; %Calculation of corrections to variables
%Balance node (node with slack number) with e and f corrections of 0
if slack==1
delt=[0;0;delt];
elseif slack==n
delt=[delt;0;0];
else
delt1=delt(1:2*slack-2,:);
delt2=delt(2*slack-1:2*n-2,:);
delt=[delt1;0;0;delt2];
end
for a=1:n
e(a,1)=e(a,1)+delt(2*a-1,1);
f(a,1)=f(a,1)+delt(2*a,1);
end
it=it+1; %Number of iterations
end
Pi=zeros(n,1); Qi=zeros(n,1);
for a=1:n %Calculating Pi, Qi
for b=1:n
Pi(a,1)=Pi(a,1)+e(a,1)*( G(a,b)*e(b,1)-B(a,b)*f(b,1) )+ f(a,1)*( G(a,b)*f(b,1)+B(a,b)*e(b,1) );
Qi(a,1)=Qi(a,1)+f(a,1)*( G(a,b)*e(b,1)-B(a,b)*f(b,1) )- e(a,1)*( G(a,b)*f(b,1)+B(a,b)*e(b,1) );
end
end
disp(['Active Power Injection at bus 2:',num2str(Pi(2,1)),' Reactive Power Injection at bus 2:',num2str(Qi(2,1))]);
figure(1);
k=1:it+1;
plot(k-1,Convergence(k,1),'-b'),xlabel('k'),ylabel('Convergence Tolerance'); hold on;
🔗 参考文献
[1]涂志军,曹晔,李伟,等.一种新型高斯塞德尔算法在电力系统潮流计算中的应用研究[J].江西科学, 2010, 28(6):4.DOI:10.3969/j.issn.1001-3679.2010.06.024.
图片
🏆团队擅长辅导定制多种科研领域MATLAB仿真，助力科研梦：
