谁能帮我指教一下改怎么编潮流计算

zhouyou529

clear;
Yb=[6.25-18.75i  -5+15i        -1.25+3.75i    0             0
    -5+15i       10.834-32.5i  -1.667+5i     -1.667+5i    -2.5+7.5i
    -1.25+3.75i  -1.667+5i     12.917-38.75i -10+30i      0
8 D% n9 V1 w/ P' s; G    0            -1.667+5i     -10+30i       12.917-38.75i  -1.25+3.75i
    0            -2.5+7.5i     0             -1.25+3.75i    3.75-11.25i];
8 W) d; r" ~/ t. ?5 h- r6 fPknow=[0;0.2;-0.45;-0.4;-0.6];Qknow=[0;0.2;-0.15;-0.05;-0.1];e1=1.06;f1=0;%已知量
e=[1.06;1;1;1;1];f=[0;0;0;0;0];
k=0;%跌代次数
! e: z6 Z9 s. z! Q$ j- T1 q4 wPQX=[1;1;1;1;1;1;1;1];%修正量   
P=zeros(4,1);Q=zeros(4,1);I=zeros(4,1);YAKEBI=zeros(8,8);
8 v/ f  u# E- Q5 kwhile (max(PQX)>=0.1|min(PQX)<=-0.1)
$ ?$ R2 T, W3 j' z/ F    k=k+1;
    for n=1:4
        for m=1:5
            P(n)=P(n)+e(n+1)*(real(Yb(n+1,m))*e(m)-imag(Yb(n+1,m))*f(m))+f(n+1)*(real(Yb(n+1,m))*f(m)+imag(Yb(n+1,m)*e(m)));
            Q(n)=Q(n)+f(n+1)*(real(Yb(n+1,m))*e(m)-imag(Yb(n+1,m))*f(m))-e(n+1)*(real(Yb(n+1,m))*f(m)+imag(Yb(n+1,m)*e(m)));
. }- F$ b- g8 F4 s        end
        PQX(2*n-1)=Pknow(n+1)-P(n);
, L9 k9 I- G; e2 A# |7 C        PQX(2*n)=Qknow(n+1)-Q(n);
$ o. I3 \7 k9 [$ Q: N7 A7 @# E        I(n)=(P(n)-Q(n)*i)/(e(n+1)-f(n+1)*i);
' x' \4 C0 `2 A7 y, C% s3 }    end
    for n=1:4
6 B4 U' w# N" u; A1 U3 y/ p        for m=1:4
            YAKEBI(2*n-1,2*m-1)=-imag(Yb(n+1,m+1))*e(n+1)+real(Yb(n+1,m+1))*f(n+1);%Hnm
: `) R& y: z/ V. r/ k            YAKEBI(2*n-1,2*m)=real(Yb(n+1,n+1))*e(n+1)+imag(Yb(n+1,n+1))*f(n+1);%Nnm
            YAKEBI(2*n,2*m-1)=-real(Yb(n+1,m+1))*e(n+1)-imag(Yb(n+1,m+1))*f(n+1);%Jnm
            YAKEBI(2*n,2*m)=-imag(Yb(n+1,m+1))*e(n+1)+real(Yb(n+1,m+1))*f(n+1);%Lnm
& r( x4 K' Z  }9 _2 L9 K. g        end
/ K) f5 o6 C+ p2 A        YAKEBI(2*n-1,2*n-1)=-imag(Yb(n+1,n+1))*e(n+1)+real(Yb(n+1,n+1))*f(n+1)+imag(I(n));%Hnn
        YAKEBI(2*n-1,2*n)=real(Yb(n+1,n+1))*e(n+1)+imag(Yb(n+1,n+1))*f(n+1)+real(I(n));%Nnn
        YAKEBI(2*n,2*n-1)=-real(Yb(n+1,n+1))*e(n+1)-imag(Yb(n+1,n+1))*f(n+1)+real(I(n));%Jnn
        YAKEBI(2*n,2*n)=-imag(Yb(n+1,n+1))*e(n+1)+real(Yb(n+1,n+1))*f(n+1)-imag(I(n));%Lnn
+ k4 G- i' A4 c# {" b    end% build the yakebi matrix
8 [  Z: x: s/ O" |    Dfe=inv(YAKEBI)*PQX; %修正fe
    for n=1:4
2 _+ x0 N% l7 \- l$ _+ z        f(n+1)=f(n+1)+Dfe(2*n-1);
        e(n+1)=e(n+1)+Dfe(2*n);
$ x& `0 n9 r( ?    end
    f
    e
end

6 e( f+ s9 [# p" Q7 X% B$ {
/ G# B+ k& G. Z  p运行会出现
Warning: Matrix is singular to working precision.
% I# ]1 c  G" y(Type "warning off MATLAB:singularMatrix" to suppress this warning.)
> In f:\MATLAB6p5\work\yy.m at line 35
/ T9 G8 E3 w5 ^4 X  M0 Q+ z
& k: ^* i+ T5 z3 b/ Q+ c0 z& [我这个程序是按照陈珩那本稳态分析的例子编的，想试试。完全没考虑PV节点
这个程序要怎么改进啊？

wangmuyidong

我没有时间仔细给你检查，但是想告诉你，你的雅克比矩阵的形成太麻烦了，极有可能存在错误，并且你得程序也不通用，雅克比矩阵的形成要仔细分析其各个表达式的关系，可以用矩阵运算直接形成。而没有必要用循环。这样会非常耗时的。