潮流计算的matlab程序实现方法

这是一个电气狗熬两个礼拜图书馆的成果,根据华中科技大学《电力系统分析》中原理编写,可用牛顿-拉夫逊和PQ分解法计算给定标幺值条件的潮流。本人水平有限,仅供参考,欢迎一起找Bug。

2019/11/17 添加算例系统图和基础数据、参考文献。

2019/01/05 添加word文档 潮流计算课程设计

2018/07/06 说明:由于本人变压器建模与PSASP不同,本人使用模型如下图,参数输入时请按该模型计算。

2018/06/18 主程序更新:增加补偿电容参数

2018/06/18 增加下载地址,详见文章底部。

应用算例系统图

基础数据(标幺值)

(1)母线数据

母线名

1

2

3

4

5

6

7

8

9

基准电压

230

230

230

230

230

230

18

13.8

16.5

(2)交流线数据

I侧母线

J侧母线

电阻

电抗

电纳的1/2

6

1

0.01

0.085

0.088

1

4

0.032

0.161

0.153

4

3

0.0085

0.072

0.0745

3

5

0.0119

0.1008

0.1045

5

2

0.039

0.17

0.179

2

6

0.017

0.092

0.079

(3)变压器数据

I侧母线

J侧母线

电抗

变比

9

6

0.0576

1

7

4

0.0625

1

8

5

0.0586

1

(4)发电机数据

母线名

母线类型

有功功率

无功功率

电压幅值

电压相角

9

slack

-

-

1.04

0

7

PV

1.63

-

1.025

-

8

PV

0.85

-

1.025

-

(5)负荷数据

母线名

母线类型

有功功率

无功功率

1

PQ

1.25

0.5

2

PQ

0.9

0.3

3

PQ

1

0.35

主程序

% file name:chaoliu_lj.m
% auther: 山东科技大学 罗江
% function:使用牛顿-拉夫逊法、PQ分解法计算潮流
% update:2018/6/18 13:22 增加补偿电容参数
%节点类型 	标号
%PQ节点 	  1
%PV节点 	  2
%slack节点  3
%能计算给定标幺值网络,有且仅有一个平衡节点的潮流
%注意:母线标号顺序要求:PQ节点-PV节点-平衡节点
%若某元件不存在,其导纳为0,阻抗为inf


clear %清除工作空间变量
clc %清屏
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%数据输入(标幺值)
SB=100; %基准容量,单位MVA

%母线基准电压
Bus=[230 230 230 230 230 230 18 13.8 16.5];

%交流线参数:I侧母线 J侧母线 阻抗 1/2接地导纳
Line=[6 1 0.01+0.085i  0.088i;
	  1 4 0.032+0.161i 0.153i; 
	  4 3 0.0085+0.072i 0.0745i;
	  3 5 0.0119+0.1008i 0.1045i;
	  5 2 0.039+0.17i 0.179i;
	  2 6 0.017+0.092i 0.079i];

%变压器参数:I侧母线 J侧母线 阻抗 变比 %变压器阻抗归算到I侧
Trans=[9 6 0.0576i 1;
	   7 4 0.0625i 1;
	   8 5 0.0586i 1];

%加接地电容器补偿: 母线 导纳
Cap=[2 0.0i];

%发电机参数:母线 节点类型 P V/U θ
Gen=[9 3 1.04 0;
	 7 2 1.63 1.025;
	 8 2 0.85 1.025]; 

%负荷参数:母线 节点类型 P Q
%按参考方向,发电机发出功率(正值),负荷消耗功率(负值)
Load=[1 1 -1.25 -0.5;
	  2 1 -0.9 -0.3;
	  3 1 -1 -0.35];

mode=1; %1-极坐标下牛拉法, 2-PQ分解法
Tmax=10; %最大迭代次数
limit=1.0e-8; %要求精度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%变压器π型等效阻抗参数
Zt=zeros(size(Trans,1),3);
Zt(:,1)=Trans(:,3)./Trans(:,4);
Zt(:,2)=Trans(:,3)./(1-Trans(:,4));
Zt(:,3)=Trans(:,3)./(Trans(:,4).^2-Trans(:,4));
Trans_pi=[Trans(:,1:2) Zt(:,1) 1./Zt(:,2) 1./Zt(:,3)];

n=numel(Bus); %总节点数
m=n-1; %PQ节点数
for i=1:size(Gen,1)%数组行数
	if Gen(i,2)==2 %除去PV节点就是PQ节点
		m=m-1;
	end
end
for i=1:size(Load,1)
	if Load(i,2)==2
		m=m-1;
	end
end
%PQ节点包含浮游节点,其PQ=0

%提取P,Q,U向量
P=zeros(1,n); %P,Q为原始数据,Pi,Qi为计算结果
Q=zeros(1,n);
U=ones(1,n); %电压初始值由此确定
cita=zeros(1,n); %此处未知节点皆设为1.0∠0 %注意:此处角度单位为度,提取后再转换成弧度,后面计算使用弧度
for i=1:size(Gen,1)
	if Gen(i,2)==1 %PQ节点
		P(Gen(i,1))=Gen(i,3);
		Q(Gen(i,1))=Gen(i,4);
	end
	if Gen(i,2)==2 %PV节点
		P(Gen(i,1))=Gen(i,3);
		U(Gen(i,1))=Gen(i,4);
	end
	if Gen(i,2)==3 %slack节点
		U(Gen(i,1))=Gen(i,3);
		cita(Gen(i,1))=Gen(i,4);
	end	
end
for i=1:size(Load,1)
	if Load(i,2)==1 %PQ节点
		P(Load(i,1))=Load(i,3);
		Q(Load(i,1))=Load(i,4);
	end
	if Load(i,2)==2 %PV节点
		P(Load(i,1))=Load(i,3);
		U(Load(i,1))=Load(i,4);
	end
	if Load(i,2)==3 %slack节点
		U(Load(i,1))=Load(i,3);
		cita(Load(i,1))=Load(i,4);
	end	
end
disp('初始条件:')
disp('各节点有功:')
disp(P);
disp('各节点无功:')
disp(Q);
disp('各节点电压幅值:')
disp(U);
cita=(deg2rad(cita)); %角度转换成弧度
disp('各节点电压相角(度):')
disp(rad2deg(cita)); %显示依然使用角度

%节点导纳矩阵的计算
Y=zeros(n); %新建节点导纳矩阵
y=zeros(n); %网络中的真实导纳
%计算y(i,j)
for i=1:size(Line,1) %与交流线联结的真实导纳
	ii=Line(i,1); jj=Line(i,2);
	y(ii,jj)=1/Line(i,3);
	y(jj,ii)=y(ii,jj);
end
for i=1:size(Trans_pi,1) %与变压器联结的真实导纳
	ii=Trans_pi(i,1); jj=Trans_pi(i,2);
	y(ii,jj)=1/Trans_pi(i,3);
	y(jj,ii)=y(ii,jj);
end
%计算y(i,i)
for i=1:size(Line,1) %与交流线联结的对地导纳
	ii=Line(i,1); jj=Line(i,2);
	y(ii,ii)=y(ii,ii)+Line(i,4);
	y(jj,jj)=y(jj,jj)+Line(i,4);
end
for i=1:size(Trans_pi,1) %与变压器联结的对地导纳
	ii=Trans_pi(i,1); jj=Trans_pi(i,2);
	y(ii,ii)=y(ii,ii)+Trans_pi(i,4);
	y(jj,jj)=y(jj,jj)+Trans_pi(i,5);
end
%算上补偿电容
for i=1:size(Cap,1)
	ii=Cap(i,1);
	y(ii,ii)=y(ii,ii)+Cap(i,2);
end
%由y计算Y
ysum=sum(y,1); %每一行求和
for i=1:n
	for j=1:n
		if i==j
			Y(i,j)=ysum(i);
		else
			Y(i,j)=-y(i,j);
		end
	end
end
disp('节点导纳矩阵:');
disp(Y);

G=real(Y); %电导矩阵
B=imag(Y); %电纳矩阵

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%以上为基础数据整理
%接下来是牛拉法的大循环
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算功率不平衡量
[dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita );
disp('有功不平衡量:');
disp(dP);
disp('无功不平衡量:');
disp(dQ);

for i=1:Tmax
    fprintf('第%d次迭代\n',i);
	%雅可比矩阵的计算
    if(mode==1)
		J=Jacobi( n,m,U,cita,B,G,Pi,Qi );
		disp('雅可比矩阵');
		disp(J);
    end

	%求解节点电压修正量
	if(mode==1)
		[dU,dcita]=Correct( n,m,U,dP,dQ,J );   
	else
    	[dU,dcita]=PQ_LJ( n,m,dP,dQ,U,B );
	end
	disp('电压、相角修正量:');
	disp(dU);
	disp(rad2deg(dcita));

	%修正节点电压
	U=U+dU;
	cita=cita+dcita;
	disp('节点电压幅值:');
	disp(U);
	disp('节点电压相角:');
	disp(rad2deg(cita));

	%计算功率不平衡量
	[dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita );
	disp('有功不平衡量:');
	disp(dP);
	disp('无功不平衡量:');
	disp(dQ);
	if (max(abs(dP))<limit && max(abs(dQ))<limit )
		break;
	end%if
end%for

%迭代结束,判断收敛
if (max(abs(dP))<limit && max(abs(dQ))<limit )
	disp('计算收敛');
else
	disp('计算不收敛或未达到要求精度');
end
%打印功率
fprintf('迭代总次数:%d\n', i);
disp('节点电压幅值:');
disp(U);
disp('节点电压相角:');
disp(rad2deg(cita));
disp('有功计算结果:');
disp(Pi);
disp('无功计算结果:');
disp(Qi);

 

子程序一

 

% filename:Unbalanced.m
% author: 山东科技大学 罗江
% function: 计算功率不平衡量
function [ dP,dQ,Pi,Qi ] = Unbalanced( n,m,P,Q,U,G,B,cita )
%计算ΔPi有功的不平衡量
for i=1:n
	for j=1:n
		Pn(j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j)));
	end
	Pi(i)=sum(Pn);
end
dP=P(1:n-1)-Pi(1:n-1); %dP有n-1个

%计算ΔQi无功的不平衡量
for i=1:n
	for j=1:n
		Qn(j)=U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j)));
	end
	Qi(i)=sum(Qn);
end
dQ=Q(1:m)-Qi(1:m); %dQ有m个


end%func

子程序二

% filename:Jacobi.m
% author:山东科技大学 罗江
% function: 计算雅可比矩阵
function [ J ] = Jacobi( n,m,U,cita,B,G,Pi,Qi )
%雅可比矩阵的计算
%分块 H N K L
%i!=j时
for i=1:n-1
	for j=1:n-1
		H(i,j)=-U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j)));
	end
end
for i=1:n-1
	for j=1:m
		N(i,j)=-U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j)));
	end
end
for i=1:m
	for j=1:n-1
		K(i,j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j)));
	end
end
for i=1:m
	for j=1:m
		L(i,j)=-U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j)));
	end
end
%i==j时
for i=1:n-1
	H(i,i)=U(i).^2*B(i,i)+Qi(i);
end
for i=1:m
	N(i,i)=-U(i).^2*G(i,i)-Pi(i);
end
for i=1:m
	K(i,i)=U(i).^2*G(i,i)-Pi(i);
end
for i=1:m
	L(i,i)=U(i).^2*B(i,i)-Qi(i);
end

%合成雅可比矩阵
J=[H N;K L];


end

子程序三

% filename:Correct.m
% author:山东科技大学 罗江
% function:修正节点电压
function [ dU,dcita ] = Correct( n,m,U,dP,dQ,J )
%求解节点电压修正量
for i=1:m
	Ud2(i,i)=U(i);
end
dPQ=[dP dQ]';
dUcita=(-inv(J)*dPQ)';
dcita=dUcita(1:n-1);
dcita=[dcita 0];
dU=(Ud2*dUcita(n:n+m-1)')';
dU=[dU zeros(1,n-m)];


end

子程序四

% filename:PQ_LJ.m
% author:山东科技大学 罗江
% function:使用PQ分解法计算电压修正量
function [ dU,dcita ] = PQ_LJ( n,m,dP,dQ,U,B )
dP_U=dP./U(1:n-1);
dQ_U=dQ./U(1:m);
dUdcita=(-inv(B(1:n-1,1:n-1))*dP_U')';
dcita=dUdcita./U(1:n-1);
dU=(-inv(B(1:m,1:m))*dQ_U')';
dU=[dU zeros(1,n-m)];
dcita=[dcita 0];%补零


end

程序下载地址:点击下载程序

参考文献

[1]《电力系统分析(第四版)》 何仰赞 温增银,华中科技大学出版社,2017.

[2]MATLAB/Simulink电力系统建模与仿真(第2版)》 于群 曹娜,机械工业出版社,2017.

[3]《电力系统分析课程设计指导及示例分析》郭丽萍 顾秀芳,中国水利水电出版社,2016.

[4]《电力系统分析的计算机算法》邱晓燕 刘天琪 黄媛,中国电力出版社,2015.

评论 59 您还未登录,请先 登录 后发表或查看评论
©️2022 CSDN 皮肤主题:大白 设计师:CSDN官方博客 返回首页

打赏作者

北窗北川

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值