潮流程序matlab
clear %清除工作空间原有的变量clc
a=textread('IEEE14.dat'); %读数据文件
n=a(1,1); %节点数
line=a(1,2); %支路数
blance=a(1,3); %平衡节点
SB=a(1,4); %基准功率
kmax=a(1,5); %最大迭代次数
ab=a(2,1); %计算精度
t=find(a(:,1)==0);
xianlu=a(,:); %形成线路参数矩阵
jiedi=a(,:);%接地支路参数矩阵
bianya=a(,:); %形成变压器支路参数矩阵
jiedian=a(,:);%节点功率数据矩阵
pv=a(,:); %pv节点数据矩阵
%读线路参数
linei=xianlu(:,2)';
linej=xianlu(:,3)';
r=xianlu(:,4)';
x=xianlu(:,5)';
b0=xianlu(:,6)';
g=r./(r.*r+x.*x);
b=-x./(r.*r+x.*x);
G=sparse(linei,linej,-g,n,n)+sparse(linej,linei,-g,n,n)+sparse(linei,linei,g,n,n)+sparse(linej,linej,g,n,n);
B=sparse(linei,linej,-b,n,n)+sparse(linej,linei,-b,n,n)+sparse(linei,linei,b+b0,n,n)+sparse(linej,linej,b+b0,n,n);
% 读变压器参数
linei=bianya(:,2)';
linej=bianya(:,3)';
r=bianya(:,4)';
x=bianya(:,5)';
k0=bianya(:,6)';
g=r./(r.*r+x.*x);
b=-x./(r.*r+x.*x);
G=G+sparse(linei,linej,-g./k0,n,n)+sparse(linej,linei,-g./k0,n,n)+sparse(linei,linei,g./k0./k0,n,n)+sparse(linej,linej,g,n,n);
B=B+sparse(linei,linej,-b./k0,n,n)+sparse(linej,linei,-b./k0,n,n)+sparse(linei,linei,b./k0./k0,n,n)+sparse(linej,linej,b,n,n);
% 读接地支路参数
jiedidian=jiedi(:,1)';
branchib=jiedi(:,2)';
B=B+sparse(jiedidian,jiedidian,branchib,n,n);
%节点导纳矩阵
Y=sparse(G+i*B);
% 读节点功率数据
pg=jiedian(:,2)/SB;
pd=jiedian(:,4)/SB;
qg=jiedian(:,3)/SB;
qd=jiedian(:,5)/SB;
p=zeros(n,1);
q=zeros(n,1);
p(jiedian(:,1))=pg-pd;
q(jiedian(:,1))=qg-qd;
% 读pv节点数据
pvjiedian=pv(:,1);
vpv=pv(:,2);
v=ones(n,1); %电压初始值
v(pvjiedian)=vpv; %pv节点(及平衡节点)电压已知
delt=zeros(n,1); %相角初始值
deltp=zeros(n,1); % 失配有功功率初值
deltq=zeros(n,1); % 失配无功功率初值
%迭代求潮流计算
for diedai=1:10 %置最大循环次数
clear i
% 失配功率方程
V=v.*cos(delt)+i*v.*sin(delt);
a2=conj(Y*V).*V;
deltp=p-real(a2);
deltq=q-imag(a2);
% 求雅可比矩阵
for i=1:n
for j=1:n
if i~=j
H(i,j)=-v(i)*v(j)*(G(i,j)*sin(delt(i)-delt(j))-B(i,j)*cos(delt(i)-delt(j)));
else
H(i,i)=v(i)*v(i)*B(i,i)+q(i)-deltq(i);
end
end
end
for i=1:n
for j=1:n
if i~=j
N(i,j)=-v(i)*(G(i,j)*cos(delt(i)-delt(j))+B(i,j)*sin(delt(i)-delt(j)));
else
N(i,i)=-v(i)*G(i,i)-(p(i)-deltp(i))/v(i);
end
end
end
for i=1:n
for j=1:n
if i~=j
J(i,j)=v(i)*v(j)*(G(i,j)*cos(delt(i)-delt(j))+B(i,j)*sin(delt(i)-delt(j)));
else
J(i,i)=v(i)*v(i)*G(i,i)-(p(i)-deltp(i));
end
end
end
for i=1:n
for j=1:n
if i~=j
L(i,j)=-v(i)*(G(i,j)*sin(delt(i)-delt(j))-B(i,j)*cos(delt(i)-delt(j)));
else
L(i,i)=v(i)*B(i,i)-(q(i)-deltq(i))/v(i);
end
end
end
%修正雅克比矩阵
Ya=;
Ya(,:)=0;
Ya(:,)=0;
Ya(blance,blance)=1; %平衡节点
Ya(blance+n,blance+n)=1;
Ya(pvjiedian+n,:)=0;
Ya(:,pvjiedian+n)=0;
for h=1:length(pvjiedian)
Ya(pvjiedian(h)+n,pvjiedian(h)+n)=1;
end
%
deltp(blance)=0; %平衡节点deltp=0
deltq(blance)=0; %平衡节点deltq=0
deltq(pvjiedian)=0; %pv节点deltq=0
deltpq=;
if max(abs(deltpq))<ab % 如果结果满足精度
break % 退出循环
end
R=Ya\(-deltpq);
delt=delt+R();
deltt=delt*180/pi
v=v+R()
end
数据格式
14201100200.1
1.00E-061
1120.019380.059170.0264
2150.054030.223040.0246
3230.046990.197970.0219
4240.058110.176320.017
5250.056950.173880.0173
6340.067010.171030.0064
7450.013350.042110
116110.094980.19890
126120.122910.255810
136130.066150.130270
147800.176150
157900.110010
169100.031810.08450
1912130.220920.199880
2013140.170930.348020
41490.127110.270380
510110.082050.192070
0
90.19
0
14700.209120.9780.9
24900.556180.969
0.91.1
35600.252020.9320.9
0
160000
26542.421.712.7
3023.3994.219
40047.8-3.9
5007.61.6
68512.2411.27.5
70000
8017.3600
90029.516.6
100095.8
11003.51.8
12006.11.6
130013.55.8
140014.95
0
11.06-4050
21.045-4050
31.01040
61.07-3040
81.09-3045
0
11052.450.00550200
244.43.510.00520100
640.63.890.00520100
0
0
0
很有用！jixuzheyangyigechengxu a! 希望能有帮助，呵呵 修正雅可比矩阵没有看懂
页:
[1]
