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

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 ]9 O. y( K4 _% Y# @" ?    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];
) F( L4 x" }  g% N6 t" I  vPknow=[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];
7 V- I$ Z' A1 U, F, Bk=0;%跌代次数
, s, O7 x1 `1 f# r" ?  PPQX=[1;1;1;1;1;1;1;1];%修正量   
P=zeros(4,1);Q=zeros(4,1);I=zeros(4,1);YAKEBI=zeros(8,8);
while (max(PQX)>=0.1|min(PQX)<=-0.1)
    k=k+1;
8 C7 M% ^" q6 g' w    for n=1:4
        for m=1:5
0 t6 f1 R' I. a7 |( H            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)));
& ?3 v/ ?2 {7 O: C/ M! `# H  p            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)));
& ^. j) }* h& p& Z2 ?        end
" h+ q3 Q* J5 L! B% m7 q        PQX(2*n-1)=Pknow(n+1)-P(n);
5 n( c6 P2 {+ K        PQX(2*n)=Qknow(n+1)-Q(n);
        I(n)=(P(n)-Q(n)*i)/(e(n+1)-f(n+1)*i);
& |5 A( s" \+ i5 D    end
    for n=1:4
3 [5 S" Z" }% C  ?. T# {        for m=1:4
$ ^5 ?7 Q9 z. L. O            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
3 |4 b; G* l6 {! ]1 I1 S! V            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
; K2 v6 U; }0 q( t& [3 q7 `        end
        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
2 h' v3 T0 o& R( j) c        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
8 V9 A0 ^  A/ D0 z4 j        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
% F- m$ m2 j8 w& b! N$ H        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
4 I( C' J( z+ z. d7 Y; q    end% build the yakebi matrix
    Dfe=inv(YAKEBI)*PQX; %修正fe
. W4 a% ~) G3 I$ N6 T    for n=1:4
        f(n+1)=f(n+1)+Dfe(2*n-1);
6 B1 u6 Q* _' `# \- S+ N        e(n+1)=e(n+1)+Dfe(2*n);
* h/ @* u7 R+ n% K: q3 j    end
    f
2 O1 I# s3 x' R    e
end
9 H$ P' r8 I/ G3 z8 b

: l: I0 ^- H1 k) b- ^运行会出现
Warning: Matrix is singular to working precision.
(Type "warning off MATLAB:singularMatrix" to suppress this warning.)
> In f:\MATLAB6p5\work\yy.m at line 35
, Z) g$ a) B' r6 q6 V6 y
我这个程序是按照陈珩那本稳态分析的例子编的，想试试。完全没考虑PV节点
& t6 z6 p3 }' b; u: \5 k- y& X. M这个程序要怎么改进啊？

wangmuyidong

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