请帮忙看看我的matlab程序错哪儿了？潮流计算

朗朗 · 发表于 2009-1-12 13:35:47

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

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

×

我对程序严重的感冒

下面是我花了2个通宵修改的程序，可是还是有问题，但是我看不出来~~请高手帮忙解决一下，谢谢
（我直接粘贴不知道行不行，如果不行后面有附件）我的邮箱langlang00000@sina.com
n=5;
nl=4;
swing=1;
pr=1e-6;
B1=[1 2 0.12i    0 1 0;%[首节点号末节点号支路阻抗对地导纳(b) 变比（无变压器则为1）是否有变压器（是为1，否为0]
2 3 0.01+0.12i 0.04 1 0;
2 4 0.2i    0 1.05 1;
2 5 0.12i    0 1 0]
B2=[1 1.02 0  0 0  0 0 0; %[节点号电压幅值电压相角发电机有功发电机无功负荷有功负荷无功节点类型（平衡节点为0，PV节点为1，PQ节点为2）]
2 0 0  0 0  0 0 2;
3 0 0  0 0 0.9 0.6  2;
4 1.0  0  1.0 0  0 0 1;
5 0 0  0 0 0.8 0.5  2]
X=[1 0;%[节点号导纳]
2 0;
3 0.1i;
4 0;
5 0.1i]
Y=zeros(n);  %初始化节点导纳矩阵
for i=1:n
if X(i,2)~=0;
      p=X(i,1);
      Y(p,p)=X(i,2);  %写入节点对地导纳
end
end
for i=1:nl
if B1(i,6)==0
      p=B1(i,1);q=B1(i,2);
else p=B1(i,2);q=B1(i,1);  %确定变压器首末节点
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));
Y(q,p)=Y(p,q);                                  %互导纳
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;  %首节点自导纳
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;             %末节点自导纳
end
Y
MAG=abs(Y);PH=angle(Y);  %求节点导纳矩阵各元素的幅值和相角

%--------------------------电压及功率初始化模块------------------------------

delta=zeros(1,n);U=zeros(1,n);
for i=1:n
delta(i)=B2(i,3);  %节点电压相角初始化
U(i)=B2(i,2);    %节点电压幅值初始化
end

P=zeros(1,n);Q=zeros(1,n);
for i=1:n
if B2(i,8)~=0
   P(i)=B2(i,4)-B2(i,6);  %节点注入有功功率初始化
end
if B2(i,8)==2
   Q(i)=B2(i,5)-B2(i,7);  %节点注入无功功率初始化
end
end

%------------------------求各节点功率不平衡量模块----------------------------

NUM=0;IT2=1;  %定义循环次数，循环条件标志
while IT2~=0
IT2=0;t1=1;t2=1;
for i=1:n

         C(i)=0;
         D(i)=0;
         for j=1:n
            C(i)=C(i)+U(i)*(U(j)*MAG(i,j)*cos(PH(i,j)+delta(j)-delta(i)));  %各节点有功功率
            D(i)=D(i)-U(i)*(U(j)*MAG(i,j)*sin(PH(i,j)+delta(j)-delta(i)));  %各节点无功功率
         end
      if  i~=swing
         DP(t1)=P(i)-C(i);    %PV节点和PQ节点的有功功率失配量
         t1=t1+1;
         if  B2(i,8)==2
            DQ(t2)=Q(i)-D(i);  %PQ节点的无功功率失配量
            t2=t2+1;
         end
      end
end

t1=t1-1;t2=t2-1;
DPQ=[DP';DQ'];  %功率失配量矩阵
for i=1:t1+t2
      if abs(DPQ(i))>pr  %收敛精度判定

         IT2=IT2+1;    %不符合精度要求，进入下一次迭代
      end
end

%---------------------------求分块雅可比矩阵模块-----------------------------

H=zeros(n);
N=zeros(n);
K=zeros(n);
L=zeros(n);  %初始化分块矩阵
for i=1:n
for j=1:n
   if i==j
         H(i,i)=-D(i)-U(i)^2*MAG(i,i)*sin(PH(i,i));
         N(i,i)= C(i)+U(i)^2*MAG(i,i)*cos(PH(i,i));
         K(i,i)= C(i)-U(i)^2*MAG(i,i)*cos(PH(i,i));
         L(i,i)= D(i)-U(i)^2*MAG(i,i)*sin(PH(i,i));                %各n阶分块矩阵对角元
   else
      H(i,j)=-U(i)*U(j)*MAG(i,j)*sin(PH(i,j)+delta(j)-delta(i));
      N(i,j)= U(i)*U(j)*MAG(i,j)*cos(PH(i,j)+delta(j)-delta(i));
      K(i,j)=-U(i)*U(j)*MAG(i,j)*cos(PH(i,j)+delta(j)-delta(i));
      L(i,j)=-U(i)*U(j)*MAG(i,j)*sin(PH(i,j)+delta(j)-delta(i));  %各n阶矩阵非对角元
   end
end
end

%----------------------------求雅可比矩阵模块-------------------------------

J=zeros(2*n);  %初始化雅可比矩阵
for i=1:n
for j=1:n
      J(i,j)=H(i,j);
      J(i,(j+n))=N(i,j);
      J((i+n),j)=K(i,j);
      J((i+n),(j+n))=L(i,j);  %将各个分块矩阵合并为2n阶雅可比矩阵
end
end

PV=[];
for i=1:n
if  B2(i,8)==1
      PV=[PV B2(i,1)];  %记录PV节点的标号
end
end

J([PV+n,swing,swing+n],:)=[];  %删除与平衡节点对应的两行，与PV节点对应的一行
J(:,[PV+n,swing,swing+n])=[];  %删除与平衡节点对应的两列，与PV节点对应的一列

J;  %最终的雅可比矩阵

%------------------------解修正方程求各节点电压模块--------------------------

modify=inv(J)*DPQ;       %各变量的修正量
Ddelta=modify([1:t1],:); %节点电压相角修正量
DU=modify([t1+1:t1+t2],:);  %节点电压幅值修正量

UR(:,NUM+1)=U(1,:);       %记录各次迭代节点电压值
t4=1;
for i=1:n
      if B2(i,8)~=0
      delta(1,i)=delta(1,i)+Ddelta(t4,1);  %修正后的节点电压相角
      t4=t4+1;
      end
end

t5=1;
for i=1:n
      if B2(i,8)==2
      U(1,i)=U(1,i)+DU(t5,1)*U(1,i);    %修正后的节点电压幅值
      t5=t5+1;
      end
end
NUM=NUM+1;  %迭代次数
if NUM==1 %最大迭代次数判断
      break;  %超过最大迭代次数，跳出
end
end

%--------------------------------输出模块----------------------------------

disp('------------------------------------------------------------------');
disp('各节点电压U幅值为（节点从小到大排列）:');
disp(U);       %输出节点电压幅值
disp('------------------------------------------------------------------');
disp('各节点电压相角为（节点从小到大排列）:');
disp(delta); %输出节点电压相角
disp('------------------------------------------------------------------');
disp('迭代次数:');
disp(NUM);    %输出迭代次数

朗朗 · 发表于 2009-1-12 14:06:13

谁帮帮忙啊，急用啊··

seeout · 发表于 2009-1-13 17:53:17

这个程序好象是一本电力系统分析教材上的，输入对了应该就没错了吧

yaochilan · 发表于 2009-1-14 10:20:27

哈哈，这个东西要研究研究了。。。。。。

norika · 发表于 2009-1-14 10:21:52

原帖由 seeout 于 2009-1-13 17:53 发表

9 f$ ?$ M2 ~0 F# m8 f5 r8 A7 d! J& p( F
登录/注册后可看大图

8 ^4 r" T% A. Y) a这个程序好象是一本电力系统分析教材上的，输入对了应该就没错了吧

那本教材？介绍一下撒，技术员同志

seeout · 发表于 2009-1-14 14:02:34

孟祥萍编著电力系统分析高等教育出版社
不过说实话，个人觉得这本书除了用上matlab有些特色外，显得粗糙了些

norika · 发表于 2009-1-14 14:09:32

原帖由 seeout 于 2009-1-14 14:02 发表
. g; H% A0 Y3 ~& N% \, J0 n* @

登录/注册后可看大图
1 q9 O: t, s5 D/ N9 L1 E" j( ?
孟祥萍编著电力系统分析高等教育出版社
5 I ]5 J/ l9 M- l9 L! y不过说实话，个人觉得这本书除了用上matlab有些特色外，显得粗糙了些

哦！我知道这本书，呵呵，不过没怎么看，谢谢你的提醒

电力新兵 · 发表于 2009-1-15 20:33:25

大家看过的书可真不少啊，俺没看过，谢谢了啊。

damayi · 发表于 2009-1-15 23:56:42

原帖由朗朗于 2009-1-12 13:35 发表

' A! V" D% u: M& V2 |& Y2 i
登录/注册后可看大图

4 T. E" _ S. n; Y/ B我对程序严重的感冒下面是我花了2个通宵修改的程序，可是还是有问题，但是我看不出来~~请高手帮忙解决一下，谢谢
, m u4 I% w u. A5 c（我直接粘贴不知道行不行，如果不行后面有附件）我的邮箱langlang00000@sina.com8 b& e, M) o# B) E. j$ b. P
n=5;4 g1 a* Y. Y7 A6 X- S& O
...

这个matlab程序写得太初级了，白白浪费了matlab的强大矩阵计算功能。简直成了C语言的样子。

^_^ · 发表于 2009-1-19 09:36:26

原帖由 damayi 于 2009-1-15 23:56 发表

* a8 j! [1 b7 p4 S$ a
登录/注册后可看大图

; T- q( ~3 A4 M$ i! A# p) [) @( V' u
2 a4 s A$ C; |- Y9 g1 B; a这个matlab程序写得太初级了，白白浪费了matlab的强大矩阵计算功能。简直成了C语言的样子。

有些这样的味道！~！

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

请帮忙看看我的matlab程序错哪儿了？潮流计算

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

相关帖子