电力系统潮流（Matlab实现）

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

5.2 负荷增加

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

5.3 发电机出力增加

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

6 MATLAB编程实例

---

6.1 潮流计算程序流程图

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

6.2 潮流计算代码

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)+j1/2line(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 (K0)&(J0) % 对地支路: 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)+KKYt;

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 (K0)&(J0) % 对地支路: 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))<EPS)&(max(abs(dAng))<EPS)

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

自我介绍一下，小编13年上海交大毕业，曾经在小公司待过，也去过华为、OPPO等大厂，18年进入阿里一直到现在。

深知大多数Python工程师，想要提升技能，往往是自己摸索成长或者是报班学习，但对于培训机构动则几千的学费，着实压力不小。自己不成体系的自学效果低效又漫长，而且极易碰到天花板技术停滞不前！

因此收集整理了一份《2024年Python开发全套学习资料》，初衷也很简单，就是希望能够帮助到想自学提升又不知道该从何学起的朋友，同时减轻大家的负担。

既有适合小白学习的零基础资料，也有适合3年以上经验的小伙伴深入学习提升的进阶课程，基本涵盖了95%以上前端开发知识点，真正体系化！

由于文件比较大，这里只是将部分目录大纲截图出来，每个节点里面都包含大厂面经、学习笔记、源码讲义、实战项目、讲解视频，并且后续会持续更新

如果你觉得这些内容对你有帮助，可以扫码获取！！！（备注Python）

的负担。**

[外链图片转存中…(img-BiuiteCO-1712799355665)]

[外链图片转存中…(img-F9ldVtec-1712799355665)]

既有适合小白学习的零基础资料，也有适合3年以上经验的小伙伴深入学习提升的进阶课程，基本涵盖了95%以上前端开发知识点，真正体系化！

由于文件比较大，这里只是将部分目录大纲截图出来，每个节点里面都包含大厂面经、学习笔记、源码讲义、实战项目、讲解视频，并且后续会持续更新

如果你觉得这些内容对你有帮助，可以扫码获取！！！（备注Python）