% file name:chaoliu14.m
% function:使用牛顿-拉夫逊法、PQ分解法计算潮流
%节点类型 标号
%PQ节点 1
%PV节点 2
%slack节点 3
%能计算给定标幺值网络，有且仅有一个平衡节点的潮流
%注意：母线标号顺序要求：PQ节点-PV节点-平衡节点
%若某元件不存在，其导纳为0，阻抗为inf
clear %清除工作空间变量
clc %清屏
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%数据输入（标幺值）
SB=100; %基准容量，单位MVA
%母线基准电压
Bus=[1 1 1 1 1 1 1 1 1 1 1 1 1 1];
%交流线参数：I侧母线 J侧母线 阻抗 1/2接地导纳
Line=[1 2 0.01938+0.05917i 0.0264;
2 3 0.04699+0.19797i 0.0219;
2 4 0.05811+0.17632i 0.0187;
1 5 0.05403+0.22304i 0.0246;
2 5 0.05695+0.17388i 0.017;
3 4 0.06701+0.17103i 0.0173;
4 5 0.01335+0.04211i 0.0064;
7 8 0+0.17615i 0;
7 9 0+0.11001i 0;
9 10 0.03181+0.0845i 0;
6 11 0.09498+0.1989i 0;
6 12 0.12291+0.15581i 0;
6 13 0.06615+0.13027i 0;
9 14 0.12711+0.27038i 0;
10 11 0.08205+0.19207i 0;
12 13 0.22092+0.19988i 0;
13 14 0.17038+0.34802i 0];
%变压器参数：I侧母线 J侧母线 阻抗 变比 %变压器阻抗归算到I侧
Trans=[5 6 0.25202i 0.932;
4 7 0.20912i 0.978;
4 9 0.55618i 0.969];
%加接地电容器补偿: 母线 导纳
Cap=[9 0.19i];
%发电机参数：母线 节点类型 P V/U θ
Gen=[1 3 1.06 0;
2 2 0.4 1.045;
3 2 -0.942 1.01;
6 2 -0.112 1.07;
8 2 0 1.09];
%负荷参数：母线 节点类型 P Q
%按参考方向，发电机发出功率（正值），负荷消耗功率（负值）
Load=[2 1 0.217 0.127;
3 1 0.942 0.19;
4 1 0.478 -0.039;
5 1 0.076 0.016;
6 1 0.112 0.075;
7 1 0 0;
8 1 0 0;
9 1 0.295 0.166;
10 1 0.09 0.058;
11 1 0.035 0.018;
12 1 0.061 0.016;
13 1 0.135 0.058;
14 1 0.149 0.05];
mode=1; %1-极坐标下牛拉法， 2-PQ分解法
Tmax=50; %最大迭代次数
limit=1.0e-5; %要求精度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%变压器π型等效阻抗参数
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);
- 1
- 2
前往页