电力系统潮流计算与PowerWorld仿真（牛顿拉夫逊法和高斯赛德尔法）（Matlab实现）

💥💥💥💞💞💞欢迎来到本博客❤️❤️❤️💥💥💥
🎉作者研究：🏅🏅🏅本科计算机专业，研究生电气学硕。主要研究方向是电力系统和智能算法、机器学习和深度学习。目前熟悉python网页爬虫、机器学习、群智能算法、深度学习的相关内容。希望将计算机和电网有效结合！⭐️⭐️⭐️

🎉🎉欢迎您的到来🎉🎉

⛅⛅⛅ 📃个人主页：科研室🌈🌈🌈

📚📚📚📋所有代码目录：电气工程科研社👨‍💻👨‍💻👨‍💻

【现在公众号名字改为：荔枝科研社】

📋📋📋本文目录如下：⛳️⛳️⛳️

目录

1 概述

2 主要任务

3 主要内容

4 案例分析

5 PowerWorld仿真

5.1 正常工作

5.2 负荷增加

5.3 发电机出力增加

6 MATLAB编程实例

6.1 潮流计算程序流程图

6.2 潮流计算代码

1 概述

最初，电力系统潮流计算是通过人工手算的。后来为了适应电力系统日益发展的需要，计算机网络已经形成，为了电力系统的潮流计算提供了物质基础。电力系统潮流计算是电力系统分析计算中最基本的内容，也是的电力系统运行及设计中必不可少的工具。根据系统给定的运行条件、网络接线及元件参数，通过潮流计算可以确定各母线电压的幅值及相角、各元件中流过的功率、整个系统的功率损耗等。

本文通过介绍基于牛顿拉夫逊法和高斯赛德尔法的潮流计算，在MATLAB中对牛顿拉夫逊法的算法进行了验证，并且用PowerWorld搭建了简单的电力系统模型，对MATLAB结果加以验证，更加形象地了解实际电力系统中潮流的分布情况。

2 主要任务

（1）在电气工程领域，潮流计算、短路计算、稳定计算俗称电力系统的三大计算。

（2）高压输电网潮流的计算机算法程序设计（PQ分解法、牛顿-拉夫逊法）

（3）或中压配电网潮流的计算机算法程序设计（前推后代法、同伦延拓法等）

（4）或电力系统短路故障的计算机算法程序设计（要求不限）

3 主要内容

（1）根据电力系统网络推导电力网络数学模型，写出节点导纳矩阵；

（2）赋予各节点电压变量（直角坐标系形式）初值后，求解不平衡量；

（3）形成雅可比矩阵；

（4）求解修正量后，重新修改初值，从2开始重新循环计算；

（5）求解的电压变量达到所要求的精度时，再计算各支路功率分布、功率损耗和平衡节点功率;

（6）上机编程调试；连调；

（7）计算分析给定系统潮流分析并与手工计算结果作比较分析。

4 案例分析

如图2-1所示，该系统由两台发电机、两台变压器、三条交流输电线路以及三个负荷组成的一个具有5节点的环网，其中节点1、2、3均为PQ节点，节点4为PV节点，节点5为平衡节点。图中参数均用标幺值表示，母线1、2、3基准电压为230kV，母线4、5基准电压为13.8kV，基准功率为100MVA。

5 PowerWorld仿真

5.1 正常工作

搭建一个基于PowerWorldd的复杂模型，如图所示，同时也表征着其正常运行状态下的潮流分布情况。

5.2 负荷增加

当节点Five上的负荷增加至530MW时，系统的潮流发生较大改变，如图所示，多条线路处于过负荷运行状态下，电压水平也降低了很多，说明在过负荷下会严重影响系统的电压水平。

5.3 发电机出力增加

当节点Six上的发电机有功出力增加至800MW时，系统的潮流发生较大改变，如图所示，多条线路处于过负荷运行状态下，电压水平也稍有降低。

6 MATLAB编程实例

6.1 潮流计算程序流程图

基于牛顿-拉夫逊法的潮流计算程序流程图如图所示。

6.2 潮流计算代码

本文仅展现部分代码，全部代码见：潮流仿真（Matlab and PowerWord) (mianbaoduo.com)


function [node_result,s_result] = PowerSystem           % 潮流计算主程序  
%% 
[node] = OpenNode; 
[nn,mn] = size(node);                        % 打开数据文件.并返回node
%% 
[line] = OpenLine;
[nl,ml] = size(line);                                    % 打开数据文件.并返回line
%% 
[node,line,nPQ,nPV,nodenum,PH,PV,PQ] = Num(node,line);    % 对节点重新排序
%% 
Y = sparse(Yij(node,line))                             % 计算节点导纳矩阵
%% 
[U] = abs(Gauss_Seidel(Y,node,nPQ,nPV))            % 返回GS算法的结果，作为初值
%% 
[node_result,s_result] =Newton_Raphson(U,Y,node,nPQ,nPV,line,nodenum);     % 用牛顿拉夫逊法计算潮流结果
%% 
Result_Write(node_result,s_result,node,line);              % 把结果写入文件
 
function [node] = OpenNode
[dfile,pathname]=uigetfile('*.m','Select Node File');   % 数据文件类型为m文件，窗口标题为“Select Node File”
if pathname == 0  
    error(' you must select a valid data file')    % 如果没有选择有效文件，则出现错误提示
else  
    lfile =length(dfile);                               % 去除文件后缀  
    eval(dfile(1:lfile-2));                             % 打开文件
end
 
 
node = ...
[
       1   1.00  0.00   -1.60  -0.80  1 ;
       2    1.00  0.00   -2.00  -1.00  1 ;
       3   1.00  0.00   -3.70  -1.30  1 ;
       4    1.05  0.00    5.00   0.00  2 ;
       5    1.05  0.00    0.00   0.00  3 ];
 
function [line] = OpenLine
[dfile,pathname]=uigetfile('*.m','Select Line File');   % 数据文件类型为m文件，窗口标题为“Select Line File”
if pathname == 0  
    error(' you must select a valid data file')     % 如果没有选择有效文件，则出现错误提示
else  
    lfile =length(dfile);                               % 去除文件后缀  
    eval(dfile(1:lfile-2));                             % 打开文件
end
 
line = ...
[
        1  2   0.04    0.25    0.0    0.25    0 ;
        1  3  0.10    0.35    0.0    0.0     0 ;
        2  3  0.08    0.30    0.0    0.25    0 ;
        5  3  0.00    0.03    0.0    0.0     1.05 ;
        4  2   0.00    0.015   0.0    0.0     1.05 ];
 
function [node,line,nPQ,nPV,nodenum,PH,PV,PQ] = Num(node,line)
%% 
[nn,mn]=size(node);             % 获取nb和nl
[nl,ml]=size(line);
%% 
nPH = 0;                        % nPH为平衡节点个数  
nPV = 0;                        % nPV为PV节点个数  
nPQ = 0;                        % nPQ为PQ节点个数  
%% 
for i = 1:nn,                   % nb为总节点数    
    type= node(i,6);            % 判断节点类型
    if type == 3,       
        nPH = nPH + 1;          % 记录个数
        PH(nPH,:)=node(i,:);    % 矩阵PH计算并储存平衡节点    
    elseif type == 2,        
        nPV = nPV +1;       
        PV(nPV,:)=node(i,:);    % 矩阵PV计算并储存PV节点    
    else        
        nPQ = nPQ + 1;      
        PQ(nPQ,:)=node(i,:);    % 矩阵PQ计算并储存PQ节点    
    end  
end
%% 
node=[PQ;PV;PH];                % 对node矩阵按PQ、PV、平衡节点的顺序重新排序
%%                          
for i=1:nl      
    for j=1:2          
        for k=1:nn              
            if  line(i,j)==nodenum(k,2)
                line(i,j)=nodenum(k,1);                 
                break              
            end          
        end      
    end  
end                             % 按排序以后的节点顺序对line矩阵重新编号
 
function Y = Yij(node,line)
%% 
[nn,mn]=size(node);
[nl,ml]=size(line);
%% 
Y=zeros(nn,nn);                        % 对导纳矩阵赋初值0 
%% 
for k=1:nl    
    I=line(k,1);   
    J=line(k,2);    
    Zt=line(k,3)+j*line(k,4);          % 读入线路参数 
    if Zt~=0
        Yt=1/Zt;                       % 接地支路不计算Yt
    end
    Ym=line(k,5)+j*1/2*line(k,6);          % 计算G+B
    K=line(k,7);                       % 取变压器线路变比
    
    %% 
    if (K==0)&(J~=0)                   % 普通线路: K=0
        Y(I,I)=Y(I,I)+Yt+Ym;           % YII为自导纳
        Y(J,J)=Y(J,J)+Yt+Ym;           
        Y(I,J)=Y(I,J)-Yt;              % 互导纳
        Y(J,I)=Y(I,J);    
    end                                %求解导纳
    %% 
    if (K==0)&(J==0)                   % 对地支路: K=0,J=0,R=X=0
        Y(I,I)=Y(I,I)+Ym;              % 直接加到自导纳上
    end        
    %% 
    if K>0                             % 变压器线路:折算到i侧,K在j侧，参考书上P108        
        Y(I,I)=Y(I,I)+Yt+Ym;        
        Y(J,J)=Y(J,J)+Yt/K/K;        
        Y(I,J)=Y(I,J)-Yt/K;        
        Y(J,I)=Y(I,J);    
    end        
    %% 
    if K<0                             % 变压器线路:折算到j侧,K在i侧        
        Y(I,I)=Y(I,I)+Yt+Ym;        
        Y(J,J)=Y(J,J)+K*K*Yt;        
        Y(I,J)=Y(I,J)+K*Yt;        
        Y(J,I)=Y(I,J);    
    end
end                                  
 
function YtYm = YtYm_NR(line)         % 计算线路的等效Yt和Ym
[nl,ml]=size(line);
%% 
YtYm = zeros(nl,5);                   % 对YtYm矩阵赋初值0，和线路行数相同，共5列
YtYm(:,1:2) = line(:,1:2);            % 矩阵前两列为节点编号，后三列分别为线路等效Yt，i侧的等效Ym，j侧的等效Ym
j = sqrt(-1);                         % 用来代表i
%% 
for k=1:nl    
    I=line(k,1);     
    J=line(k,2);    
    Zt=line(k,3)+j*line(k,4);         % 读入线路参数
    if Zt~=0
        Yt=1/Zt;                      % 接地支路不计算Yt
    end
    Ym=line(k,5)+j*line(k,6);         % 计算G+B
    K=line(k,7);                      % 取变压器线路变比
    
    %% 
    if (K==0)&(J~=0)                 % 普通线路: K=0
        YtYm(k,3) = Yt;               
        YtYm(k,4) = Ym;
        YtYm(k,5) = Ym;
    end    
    %% 
    if (K==0)&(J==0)                 % 对地支路: K=0,J=0,R=X=0
        YtYm(k,4) = Ym;    
    end        
    %% 
    if K>0                           % 变压器线路: 参数折算到i侧,K在j侧，参考书上P108        
        YtYm(k,3) = Yt/K;        
        YtYm(k,4) = Ym+Yt*(K-1)/K;        
        YtYm(k,5) = Yt*(1-K)/K/K;    
    end 
    %% 
    if K<0                           % 变压器线路: 参数折算到j侧,K在i侧        
        YtYm(k,3) = -Yt*K;        
        YtYm(k,4) = Ym+Yt*(1+K);        
        YtYm(k,5) = Yt*(K^2+K);     
    end
end
 
 
%% 
while k                                                % k不为0就继续计算
    for i=1:nn
        switch node(i,6)                               % 判断节点类型 
            case 1                                     % 1为PQ节点
               [U1] = PQ(i,U,node,Y,nPQ);              % 代入PQ子程序中
               U(i,1)=U1(i,1);                         % 把值放到U矩阵中
               
            case 2                                     % 2为PV节点
               [U2] = PV(i,U,node,Y,nPV);              % 代入PV子程序中
               U(i,1)=U2(i,1);
               
            case 3                                     % 3为平衡节点，不变
               U(i,1)=node(i,2)*cos(node(i,3))-node(i,2)*sin(node(i,3))*sqrt(-1);
        end  
    end
        k=k-1;
end
 
function [node_result,s_result] =Newton_Raphson(U,Y,node,nPQ,nPV,line,nodenum)    
% 牛顿拉夫逊法主程序
[nn,mn] = size(node);
%% 
EPS = 1.0e-10;                                                 % 设定误差精度
%% 
for t = 1:1000                   % 开始迭代，最大迭代次数为1000，不收敛时跳出
%% 
[dP,dQ] = DPQ(Y,node,nPQ,nPV);                                                    % 计算P与Q的偏差
%% 
Jac = Jac_NR(node,Y,nPQ); 
%Jac = sparse(Jac_NR(node,Y,nPQ));                              % 计算雅克比矩阵
%% 
UD = zeros(nPQ,nPQ); 
for i = 1:nPQ
    UD(i,i) = U(i,1);                   % 生成电压对角矩阵，方便计算，在书上P117
end
%% 
dAngU = Jac \ [dP;dQ];                                         % 公式在P117
dAng = dAngU(1:nn-1,1);                                       % 计算相角修正量
dU = UD*(dAngU(nn:nn+nPQ-1,1));                             % 计算电压修正量
node(1:nPQ,2) = node(1:nPQ,2) - dU;                                                          % 修正电压
node(1:nn-1,3) = node(1:nn-1,3) - dAng;                            % 修正相角
%% 
if (max(abs(dU))
   break
end                                                                            % 判断是否满足精度，满足就跳出
end
%% 
node = PQ_NR(node,Y,nPQ,nPV)                     % 计算每个节点的有功和无功注入
%% 
[node,line] = ReNum(node,line,nodenum);               % 对节点恢复以前编号
%% 
YtYm = YtYm_NR(line);                  % 计算线路的等效Yt和Ym，以计算线路潮流
%% 
node_result = Node_result(node);               % 计算节点数据结果
s_result = S_result(node,line,YtYm);                    % 计算线路潮流
 
 
function U1 = PQ(i,U,node,Y,nPQ)                    % PQ节点的电压计算
%% 
[nn,mn]=size(node);                                % 获取参数和初始化矩阵
N1=zeros(nPQ,1);
N2=zeros(nn,1);
U1=zeros(nPQ,1);
%% 
N1(i,1)=((node(i,4)-node(i,5)*sqrt(-1))/(conj(U(i,1))))/Y(i,i);             % 算式前一项
%% 
for j=1:nn                                % 算式求和项
     if j~=i
        N2(i,1)=N2(i,1)+Y(i,j)*U(j,1);
     end
end
           U1(i,1)=(N1(i,1)-N2(i,1)/Y(i,i));                  % 计算最终值
 
 
function U2 = PV(i,U,node,Y,nPV)                        % 计算PV节点电压的子程序
%% 
[nn,mn]=size(node);                  % 只保存i的矩阵，只需要相应类型的节点数
N3=zeros(nn,1);
N4=zeros(nPV,1);
N5=zeros(nPV,1);
N6=zeros(nn,1);
U2=zeros(nPV,1);
%% 
for j=1:nn                                                         % 计算无功
     N3(i,1)=N3(i,1)+conj(Y(i,j))*conj(node(i,2)*cos(node(i,3))+node(i,2)*sin(node(i,3))*sqrt(-1));
end
N4(i,1)=imag(U(i,1)*N3(i,1));
N5(i,1)=((node(i,4)-N4(i,1)*sqrt(-1))/(conj(U(i,1))))/Y(i,i);
%% 
for j=1:nn                                                 % 用无功计算相角
    if j~=i
       N6(i,1)=N6(i,1)+Y(i,j)*(node(i,2)*cos(node(i,3))+node(i,2)*sin(node(i,3))*sqrt(-1));
    end
end                                                        % 保持原幅值
           U2(i,1)=node(i,2)*cos(angle(N5(i,1)-N6(i,1)))+node(i,2)*sin(angle(N5(i,1)-N6(i,1)))*sqrt(-1);
 
 
function node = PQ_NR(node,Y,nPQ,nPV)                  % 计算功率注入
%% 
n = nPQ+nPV+1;                                         % n为总节点数
%% 
for i = nPQ+1:n-1                               % 利用公式计算PV节点的无功注入
    for j = 1:n
        node(i,5)=node(i,5)+node(i,2)*node(j,2)*(real(Y(i,j))*sin(node(i,3)-node(j,3))-imag(Y(i,j))*cos(node(i,3)-node(j,3)));
    end
end                   
%% 
for j =1:n                                 % 计算平衡节点的无功及有功注入
    node(n,4)=node(n,4)+node(n,2)*node(j,2)*(real(Y(n,j))*cos(node(n,3)-node(j,3))+imag(Y(n,j))*sin(node(n,3)-node(j,3)));
    node(n,5)=node(n,5)+node(n,2)*node(j,2)*(real(Y(n,j))*sin(node(n,3)-node(j,3))-imag(Y(n,j))*cos(node(n,3)-node(j,3)));
end
 
function [node,line] = ReNum(node,line,nodenum)         % 用来恢复编码
%% 
[nn,mn]=size(node);
[nl,ml]=size(line);
node_temp = zeros(nn,mn);                               % 暂时存放用的矩阵
k = 1;                                % node_temp矩阵用于临时存放bus矩阵的数据
%% 
for i = 1 :nn   % 利用node矩阵的首列编号重新对node矩阵排序并存入node_temp矩阵中
    for j = 1 : nn
        if node(j,1) == k
            node_temp(k,:) = node(j,:);
            k = k + 1;
        end
    end
end
%% 
node = node_temp;                                       % 重新赋值回node
%% 
for i=1:nl      
    for j=1:2          
        for k=1:nn              
            if  line(i,j)==nodenum(k,1)
                line(i,j)=nodenum(k,2);                 
                break              
            end          
        end      
    end  
end                                                     % 恢复line的编号
 
function s_result = S_result(node,line,YtYm)          % 计算功率分布，主要参考P108变压器的等值电路图
%% 
[nl,ml]=size(line);
s_result = zeros(nl,5);                                     % 储存线路潮流计算结果
s_result(:,1:2) = line(:,1:2);                                % 前两列为节点
%% 
for k=1:nl
    I = s_result(k,1);
    J = s_result(k,2);                                        % 储存节点
    %%                                     
    if (J~=0)&(I~=0)                         % 计算非接地支路的潮流，分别计算i侧的功率和j侧的功率然后求和
        s_result(k,3)=node(I,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,4)))-node(I,2)*node(J,2)*(cos(node(I,3))+j*sin(node(I,3)))*(conj(cos(node(J,3))+j*sin(node(J,3))))*conj(YtYm(k,3));
        s_result(k,4)=node(J,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,5)))-node(I,2)*node(J,2)*(conj(cos(node(I,3))+j*sin(node(I,3))))*(cos(node(J,3))+j*sin(node(J,3)))*conj(YtYm(k,3));
        s_result(k,5)=s_result(k,3) + s_result(k,4);          
    elseif J==0                                   % 利用公式计算接地支路的潮流，上面是i侧，下面是j侧
        s_result(k,3)=node(I,2)^2*conj(YtYm(k,4));
        s_result(k,5)=s_result(k,3)+s_result(k,4);
    else
        s_result(k,4)=node(J,2)^2*conj(YtYm(k,5));
        s_result(k,5)=s_result(k,3)+s_result(k,4);          
    end
end
 
function node_result = Node_result(node);             % 计算节点结果
%% 
[nn,mn]=size(node);
%% 
node_result = zeros(nn,4);                            % node_result储存计算结果
node_result(:,1:2) = node(:,1:2);                     % 这两列为节点编号和电压幅值
node_result(:,3) = node(:,3) *180 / pi;               % 转化成角度
node_result(:,4) = node(:,4) + (sqrt(-1))*node(:,5);  % 注入功率
 
 
function [] = Result_Write(node_result,s_result,node,line)
[nn,mn] = size(node);
[nl,ml] = size(line);
ace = fopen('Result.TXT','a');                           % 结果保存在Result.txt中
%% 
fprintf(ace,'电力系统分析课设第四组\n节点计算结果：\n节点      节点电压        节点相角（角度）          节点注入功率\n');
for i = 1:nn
    fprintf(ace,'%2.0f      ', node_result(i,1));
    fprintf(ace,'%10.6f      ',node_result(i,2));
    fprintf(ace,'%10.6f      ',node_result(i,3));
    fprintf(ace,'%10.6f + j %10.6f\n',real(node_result(i,4)),imag(node_result(i,4)));
end
%% 
fprintf(ace,'\n线路计算结果：\n节点I    节点J            线路功率S(I,J)             线路功率S(J,I)               线路损耗dS(I,J)\n');
for i = 1:nl
    fprintf(ace,'%2.0f      ',s_result(i,1));
    fprintf(ace,'%2.0f      ',s_result(i,2));
    fprintf(ace,'%10.6f + j %10.6f     ',real(s_result(i,3)),imag(s_result(i,3)));
    fprintf(ace,'%10.6f + j %10.6f     ',real(s_result(i,4)),imag(s_result(i,4)));
    fprintf(ace,'%10.6f + j %10.6f\n',real(s_result(i,5)),imag(s_result(i,5)));  
end
fclose(ace);

在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量的分析比较供电方案或运行方式的合理性、可靠性和经济性。此外，电力系统的潮流计算也是计算机系统动态稳定和静态稳定的基础，所以潮流计算是研究电力系统的一种很重要和基础的计算。它的发展主要围绕这样几个方面：计算方法的收敛性、可靠性；计算速度的快速性；对计算机存储容量的要求以及计算的方便、灵活等。

相关阅读:
5. 最长回文子串
 Docker学习笔记
 [附源码]Python计算机毕业设计SSM科技类产品众筹系统（程序+LW)
删除的文件夹不在回收站中如何恢复呢?
【牛客网面试必刷TOP101】二分查找/排序
 Flink报错could not be loaded due to a linkage failure
Unity Xlua热更新框架（三）：资源管理
 Vue3：创建脚手架
 Java开发者的神经网络进阶指南：深入探讨交叉熵损失函数
 人体姿态估计(人体关键点检测)2D Pose训练代码和Android源码
原文地址：https://blog.csdn.net/weixin_46039719/article/details/126260451