牛拉法潮流程序

星夜飘雪 · 发表于 2010-7-15 09:31:48

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

您需要登录才可以下载或查看，没有帐号？立即加入

x

希望对大家有所帮助哦

clc
disp('此程序为牛拉法解潮流')
nPQ=input('请输入PQ节点的个数：');
nPV=input('请输入PV节点的个数：');
n=nPQ+nPV+1;
Ps=[0;-0.5;0.2];
Qs=[0;-0.3];
Us=[1.0+j*0;1.0+j*0;1.05+j*0;1.05+j*0];
% nl nr R X Bl Br
zdata=[1 2 0 0.1880 -0.6815 0.6040;
1 3 0.1302 0.2479 0.0129 0.0129;
1 4 0.1736 0.3306 0.0172 0.0172;
3 4 0.2603 0.4959 0.0259 0.0259];
% nx B
xdata=[2 0.05];
dPQU=1;
%计算导纳矩阵
nl=zdata(:,1);
nr=zdata(:,2);
R=zdata(:,3);
X=zdata(:,4);
Bl=zdata(:,5);
Br=zdata(:,6);
nx=xdata(:,1);
Bx=xdata(:,2);
nbr=length(zdata(:,1));
nbrx=length(xdata(:,1));
Z=R+j*X;
y=ones(nbr,1)./Z;
Y=zeros(n,n);
%计算非对角元素
for ii=1:nbr
Y(nl(ii),nr(ii))= Y(nl(ii),nr(ii))-y(ii);
Y(nr(ii),nl(ii))= Y(nl(ii),nr(ii));
end
%计算对角元素
for ii=1:n
for jj=1:nbr
if nl(jj)==ii|nr(jj)==ii
Y(ii,ii)=Y(ii,ii)+y(jj);
end
end
end
for ii=1:nbr
Y(nl(ii),nl(ii))= Y(nl(ii),nl(ii))+j*Bl(ii);
Y(nr(ii),nr(ii))= Y(nr(ii),nr(ii))+j*Br(ii);
end
for ii=1:nbrx
Y(nx(ii),nx(ii))=Y(nx(ii),nx(ii))+j*Bx(ii);
end
%分离G、B
G=real(Y);
B=imag(Y);
disp('导纳矩阵：');
Y
e=real(Us);
f=imag(Us);
k=0;
while dPQU>0.00001
%求dP
dP=zeros(n-1,1);
for ii=1:n-1
t=0;
for jj=1:n
t=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));
end
dP(ii)=Ps(ii)-real(t*(e(ii)+j*f(ii)));
end
%求dQ
dQ=zeros(nPQ,1);
for ii=1:nPQ
t=0;
for jj=1:n
t=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));
end
dQ(ii)=Qs(ii)-imag(t*(e(ii)+j*f(ii)));
end
%求dU^2
dU2=zeros(nPV,1);
ii=1:nPV;
dU2(ii)=abs(Us(ii+nPQ))^2-abs(e(ii+nPQ)+j*f(ii+nPQ))^2;
dS=[dP;dQ;dU2];
dPQU=max(abs(dS));
if(dPQU>0.00001)
k=k+1
%形成雅克比行列式
Jacob=zeros(2*(n-1),2*(n-1));
%P部分
for ii=1:n-1
mid1=0;
mid2=0;
for jj=1:n
if ii~=jj&&jj<n
Jacob(ii,2*jj-1)=-(G(ii,jj)*e(ii)+B(ii,jj)*f(ii));
Jacob(ii,2*jj)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);
end
mid1=mid1+G(ii,jj)*f(jj)+B(ii,jj)*e(jj);
mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj);
end
Jacob(ii,2*ii-1)=-mid2-G(ii,ii)*e(ii)-B(ii,ii)*f(ii);
Jacob(ii,2*ii)=-mid1+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);
end
%Q部分
for ii=1:nPQ
mid1=0;
mid2=0;
for jj=1:n
if ii~=jj&&jj<n
Jacob(ii+n-1,2*jj-1)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);
Jacob(ii+n-1,2*jj)=G(ii,jj)*e(ii)+B(ii,jj)*f(ii);
end
mid1=mid1+G(ii,jj)*f(jj)+B(ii,jj)*e(jj);
mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj);
end
Jacob(ii+n-1,2*ii-1)=mid1+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);
Jacob(ii+n-1,2*ii)=-mid2+G(ii,ii)*e(ii)+B(ii,ii)*f(ii);
end
%U2部分
for ii=nPQ+1:n-1
Jacob(ii+n-1,2*ii-1)=-2*e(ii);
Jacob(ii+n-1,2*ii)=-2*f(ii);
end
dU=-inv(Jacob)*dS;
de=zeros(n-1,1);
df=zeros(n-1,1);
ii=1:n-1;
de(ii)=dU(2*ii-1);
df(ii)=dU(2*ii);
e(ii)=e(ii)+de(ii);
f(ii)=f(ii)+df(ii);
end
end%迭代结束
U=e+j*f;
%计算PV节点的Q
P=zeros(n,1);
Q=zeros(n,1);
for ii=1:nPV
t=0;
for jj=1:n
t=t+conj(Y(ii+nPQ,jj))*(e(jj)-j*f(jj));
end
Q(ii+nPQ)=imag(t*(e(ii+nPQ)+j*f(ii+nPQ)));
end
%计算平衡节点
t=0;
for jj=1:n
t=t+conj(Y(n,jj))*(e(jj)-j*f(jj));
end
P(n)=real(t*(e(n)+j*f(n)));
Q(n)=imag(t*(e(n)+j*f(n)));
ii=1:n-1;
P(ii)=Ps(ii);
ii=1:nPQ;
Q(ii)=Qs(ii);
%计算线路潮流
Sij=zeros(nbr,1);
Sji=zeros(nbr,1);
dSij=zeros(nbr,1);
for ii=1:nbr
Sij(ii)=U(nl(ii))*(conj(U(nl(ii)))*(-j*Bl(ii))+(conj(U((nl(ii))))-conj(U((nr(ii)))))*conj(y(ii)));
Sji(ii)=U(nr(ii))*(conj(U(nr(ii)))*(-j*Br(ii))+(conj(U((nr(ii))))-conj(U((nl(ii)))))*conj(y(ii)));
dSij(ii)=Sij(ii)+Sji(ii);
end
nn=[1:n]';
disp(' n e f P Q');
Display1=[nn e f P Q]
disp(' nl nr Sij Sji dSij');
Display2=[nl nr Sij Sji dSij]

复制代码

leichen · 发表于 2010-12-11 19:35:42

我正好要用，虽然不知道这个能不能行，先感谢楼主分享啦

wangwangmoon · 发表于 2010-12-24 11:19:59

不行吧，不能保证雅可比矩阵是非奇异矩阵，不能直接求逆吧

花尹儿 · 发表于 2016-3-29 16:46:35

正在学习中

云来123 · 发表于 2016-3-29 17:01:51

谢谢分享

风晨 · 发表于 2018-12-27 08:50:46

原创，永远支持~~~~~~~~~~~~

帐号		自动登录	找回密码
密码			立即加入

[分享] 牛拉法潮流程序

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

相关帖子