跪求MATLAB实现的潮流计算实例

yaxuan · 发表于 2010-5-6 14:17:24

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

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

×

本人小白

希望那个好心的提供一个用MATLAB实现的潮流计算实例，谢谢啊

tonysh · 发表于 2010-5-15 17:22:43

比较权威的，建议参考Matpower，网上搜索一下

水木羽泉 · 发表于 2010-5-15 22:24:10

就在论坛里找，好多的

百草春生 · 发表于 2010-5-17 17:35:18

楼主，有了可以发我吗？谢谢wjdf1234@163.com

沧山牧水 · 发表于 2010-7-28 12:24:23

%本程序采用牛顿-拉夫讯法对一个五节点电力系统进行潮流计算
% 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*j  0    1;
2 3 0.08+0.30*j 0.25*j 1;
2 4  0.015*j 0    1.05;
3 5  0.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.05  1.05 3;
0          1.05  0 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(:,[1:NO-2]);
   DP=J([1:NO-2],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([1:8],9)');
   MDU(Ci)=max(abs(J([1:8],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(E(p))-conj(E(q))]*conj(1./(Br(i,3)));
SBS(i)=[abs(E(p))]^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(E(q))-conj(E(p))]*conj(1./(Br(i,3)));
SBE(i)=[abs(E(q))]^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');

fulepingyi · 发表于 2010-7-28 13:55:14

matpower在潮流计算方面多不错的

smalling321 · 发表于 2010-9-2 15:43:18

matpower

是什么？请指导一下。

韩小姐 · 发表于 2011-6-29 12:08:26

回复 3# 水木羽泉

找不着啊，帮帮忙，刚刚接触这东西都不知道怎么找

chongchong · 发表于 2012-2-5 13:14:19

会编程的人都是牛人呀

lhl520lh · 发表于 2012-2-19 18:22:19

前期还是要多看看别人的程序。。。然后自己去编。。。

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

跪求MATLAB实现的潮流计算实例

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

相关帖子