跪求MATLAB实现的潮流计算实例
本人小白::cry:: 希望那个好心的提供一个用MATLAB实现的潮流计算实例，谢谢啊 比较权威的，建议参考Matpower，网上搜索一下 就在论坛里找，好多的 楼主，有了可以发我吗？谢谢wjdf1234@163.com %本程序采用牛顿-拉夫讯法对一个五节点电力系统进行潮流计算% n=input('请输入系统节点数目:n=');
% nl=input('请输入系统支路数目:nl=');
% ph=input('请输入平衡节点号:ph=');
% jd=input('请输入误差精度：jd=');
% B=input('请输入由支路参数形成的矩阵:B=');
% A=input('请输入各节点参数形成的矩阵:A=');
clear all
n=5;
nl=5;
ph=5;
jd=1E-6;
%支路参数矩阵
Br=[1 2 0.04+0.25*j 0.25j 1;
1 3 0.1+0.35*j0 1;
2 3 0.08+0.30*j 0.25*j 1;
2 40.015*j 0 1.05;
3 50.03*j 0 1.05];
%节点参数矩阵:2-PQ节点，3-PV节点，4-平衡节点
A=[-1.6-0.8*j 1 0 2;
-2-1*j 1 0 2;
-3.7-1.3*j 1 0 2;
5 1.051.05 3;
0 1.050 4];
%雅克比矩阵形成
Y=zeros(n);
e=zeros(1,n);
f=zeros(1,n);
V=zeros(1,n);
for i=1:nl
p=Br(i,1);
q=Br(i,2);
Y(p,q)=Y(p,q)-1./(Br(i,3)*Br(i,5));
Y(q,p)=Y(p,q);
Y(q,q)=Y(q,q)+1./(Br(i,3))+Br(i,4);
Y(p,p)=Y(p,p)+1./(Br(i,3)*Br(i,5)^2)+Br(i,4);
end
disp('节点导纳矩阵为：Y=');
disp(Y);
G=real(Y);
B=imag(Y);
for i=1:n
e(i)=real(A(i,2));
f(i)=imag(A(i,2));
V(i)=A(i,3);
end
for i=1:n
S(i)=A(i);
end
P=real(S);
Q=imag(S);
%雅克比矩阵求取
Ci=0;
a=1;
NO=2*n;
N=NO-1;
while a~=0
a=0;
for i=1:n
if i~=ph
C(i)=0;
D(i)=0;
for p=1:n
C(i)=C(i)+G(i,p)*e(p)-B(i,p)*f(p);
D(i)=D(i)+G(i,p)*f(p)+B(i,p)*e(p);
end
P1=e(i)*C(i)+f(i)*D(i);
Q1=f(i)*C(i)-e(i)*D(i);
V2=e(i)^2+f(i)^2;
if A(i,4)~=3
DP=P(i)-P1;
DQ=Q(i)-Q1;
for k=1:n
%非平衡节点时，非对角线元素
if k~=ph & k~=i
X1=-G(i,k)*f(i)+B(i,k)*e(i);
X2=-G(i,k)*e(i)-B(i,k)*f(i);
X3=-X2;
X4=X1;
p=2*i-1;
q=2*k-1;
J(p,q)=X1;
J(p,N)=DP;
m=p+1;
J(m,q)=X3;
J(m,N)=DQ;
q=q+1;
J(p,q)=X2;
J(m,q)=X4;
%非平衡节点时对角线元素
elseif k~=ph & k==i
X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
X2=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);
X3=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);
X4=D(i)+B(i,i)*e(i)-G(i,i)*f(i);
p=2*i-1;
q=2*k-1;
J(p,q)=X1;
J(p,N)=DP;
m=p+1;
J(m,q)=X3;
J(m,N)=DQ;
q=q+1;
J(p,q)=X2;
J(m,q)=X4;
end
end
else
DP=P(i)-P1;
DV=V(i)^2-V2;
for s=1:n
if s~=ph & s~=i
X1=-G(i,s)*f(i)+B(i,s)*e(i);
X2=-G(i,s)*e(i)-B(i,s)*f(i);
X5=0;
X6=0;
p=2*i-1;
q=2*s-1;
J(p,q)=X1;
J(p,N)=DP;
m=p+1;
J(m,q)=X5;
J(m,N)=DV;
q=q+1;
J(p,q)=X2;
J(m,q)=X6;
elseif s~=ph & s==i
X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
X2=-C(i)-G(i,s)*e(i)-B(i,s)*f(i);
X5=-2*f(i);
X6=-2*e(i);
p=2*i-1;
q=2*s-1;
J(p,q)=X1;
J(p,N)=DP;
m=p+1;
J(m,q)=X5;
J(m,N)=DV;
q=q+1;
J(p,q)=X2;
J(m,q)=X6;
end
end
end
end
end
disp('系统雅克比矩阵为：J=');
disp(J);
%对于高压电力系统，最大元素出现在对角线位置，简单计算如下
DJ=J(:,);
DP=J(,N);
DU=DJ\DP;
disp('************************************');
disp(' ')
disp('第几次的修正值DU？');
disp(Ci);
disp('***********************************');
disp(DU);
for i=1:NO-2
eps=abs(DP(i));
if eps>=jd;
a=a+1;
end
end
Ci=Ci+1;
for i=1:n-1
e(i)=e(i)-DU(2*i);
f(i)=f(i)-DU(2*i-1);
end
E=e+j*f;
disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
disp(' ');
disp('第几次迭代后的近似计算值');
disp(Ci);
disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
for i1=1:n
V(i1)=sqrt(e(i1)^2+f(i1)^2);
ANGLE(i1)=atan(f(i1)./e(i1))*(180./pi);
end
%下列输出结果均按节点号从小到大排列
disp('各个节点的电压实际值为：E=');
disp(E);
disp('各个节点电压幅值为：V=');
disp(V);
disp('各个节点的电压相交为：ANGLE=');
disp(ANGLE);
disp('各个节点的功率误差变化');
disp(J(,9)');
MDU(Ci)=max(abs(J(,9)'));
for p=1:n
Sm(p)=0;
for q=1:n
Sm(p)=Sm(p)+conj(Y(p,q))*conj(E(q));
end
S(i)=E(p)*Sm(p);
end
end
disp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
disp(' ');
disp('计算的最终功率近似值如下各项值');
disp(' ');
disp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
disp('各个节点功率为:S=');
disp(S);
disp('各条支路的首端功率Si为:');
for i=1:nl
p=Br(i,1);q=Br(i,2);
A(i)=E(p)**conj(1./(Br(i,3)));
SBS(i)=^2*conj(Br(i,4))+A(i);
end
disp(SBS);
disp('各条支路的末端功率Sj为:');
for i=1:nl
p=Br(i,1);q=Br(i,2);
A2=E(q)**conj(1./(Br(i,3)));
SBE(i)=^2*conj(Br(i,4))+A2;
end
disp(SBE);
disp('各条支路的功率损耗为:');
for i=1:nl
p=Br(i,1);q=Br(i,2);
SBL(i)=SBS(i)+SBE(i);
end
disp(SBL);
%绘制功率误差曲线，先用最小二乘法进行曲线拟合
for i=1:Ci
CSH(i)=i;
end
P=polyfit(CSH,MDU,3);
plot(CSH,MDU,'-O'); matpower在潮流计算方面多不错的 matpower
是什么？请指导一下。 回复 3# 水木羽泉
找不着啊 ，帮帮忙 ，刚刚接触这东西 都不知道怎么找 会编程的人都是牛人呀 前期还是要多 看看别人的程序。。。然后自己去编。。。
页:
[1]
2
