交流社区 › MATLAB › 请帮忙看看我的matlab程序错哪儿了？潮流计算

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

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

我对程序严重的感冒::sweat:: 下面是我花了2个通宵修改的程序，可是还是有问题，但是我看不出来~~请高手帮忙解决一下，谢谢
（我直接粘贴不知道行不行，如果不行后面有附件）我的邮箱langlang00000@sina.com
n=5;
nl=4;
swing=1;
pr=1e-6;
B1=
    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=
    2 0    00   00   0   2;
    3 0    00   0 0.9 0.62;
    4 1.001.0 00   0   1;
    5 0    00   0 0.8 0.52]
X=
    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
      ifi~=swing
            DP(t1)=P(i)-C(i);      %PV节点和PQ节点的有功功率失配量
            t1=t1+1;
            ifB2(i,8)==2
                DQ(t2)=Q(i)-D(i);%PQ节点的无功功率失配量
                t2=t2+1;
            end
      end
    end
   
    t1=t1-1;t2=t2-1;
    DPQ=;%功率失配量矩阵
    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
    ifB2(i,8)==1
      PV=;%记录PV节点的标号
    end
end

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

J;%最终的雅可比矩阵

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

    modify=inv(J)*DPQ;          %各变量的修正量
    Ddelta=modify(,:);    %节点电压相角修正量
    DU=modify(,:);%节点电压幅值修正量
   
    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

谁帮帮忙啊，急用啊··::cry::

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

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

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

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

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

原帖由 seeout 于 2009-1-13 17:53 发表 https://tech.cepsc.com/images/common/back.gif
这个程序好象是一本电力系统分析教材上的，输入对了应该就没错了吧::lol:: ::lol::

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

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

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

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

原帖由 seeout 于 2009-1-14 14:02 发表 https://tech.cepsc.com/images/common/back.gif
孟祥萍 编著 电力系统分析 高等教育出版社
不过说实话，个人觉得这本书除了用上matlab有些特色外，显得粗糙了些

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

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

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

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

原帖由 朗朗 于 2009-1-12 13:35 发表 https://tech.cepsc.com/images/common/back.gif
我对程序严重的感冒::sweat:: 下面是我花了2个通宵修改的程序，可是还是有问题，但是我看不出来~~请高手帮忙解决一下，谢谢
（我直接粘贴不知道行不行，如果不行后面有附件）我的邮箱langlang00000@sina.com
n=5;
...
这个matlab程序写得太初级了，白白浪费了matlab的强大矩阵计算功能。简直成了C语言的样子。

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

原帖由 damayi 于 2009-1-15 23:56 发表 https://tech.cepsc.com/images/common/back.gif

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

有些这样的味道！~！

查看完整版本: 请帮忙看看我的matlab程序错哪儿了？潮流计算