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

zhouyou529

clear;
7 h( N: ~6 W" t+ l7 }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
9 W# l: V+ K9 l: o0 b) Q2 q6 Z    -1.25+3.75i  -1.667+5i     12.917-38.75i -10+30i      0
' ?( \8 b; x: C8 L    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];
- Y* A9 R, f# u: xPknow=[0;0.2;-0.45;-0.4;-0.6];Qknow=[0;0.2;-0.15;-0.05;-0.1];e1=1.06;f1=0;%已知量
* G& H7 e& w5 e. `8 ke=[1.06;1;1;1;1];f=[0;0;0;0;0];
: r. S  a! p4 X+ q; ?/ P9 ck=0;%跌代次数
  M( S1 Y% s8 Z5 k9 l! X2 y9 J1 @PQX=[1;1;1;1;1;1;1;1];%修正量   
4 \" i) `6 a+ l' g4 C+ k8 n5 lP=zeros(4,1);Q=zeros(4,1);I=zeros(4,1);YAKEBI=zeros(8,8);
* C$ X0 w/ A+ D9 G4 m0 e* R3 Zwhile (max(PQX)>=0.1|min(PQX)<=-0.1)
    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)));
        end
9 n( g; j8 i) P; ]7 V        PQX(2*n-1)=Pknow(n+1)-P(n);
        PQX(2*n)=Qknow(n+1)-Q(n);
7 C  B- L( J6 R" o- X        I(n)=(P(n)-Q(n)*i)/(e(n+1)-f(n+1)*i);
    end
    for n=1:4
5 I! Z$ T7 i2 z; N, r( E        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
9 E8 J1 @" i& j8 H" u- U            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
! h9 P" G. P  ]4 ^            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
/ y9 n/ a: ]! [  b) k        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
4 x" G, G3 \9 i' B" b        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
4 z$ A3 W8 r. f        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
    end% build the yakebi matrix
    Dfe=inv(YAKEBI)*PQX; %修正fe
    for n=1:4
- C: U! k5 y; N0 X. O        f(n+1)=f(n+1)+Dfe(2*n-1);
        e(n+1)=e(n+1)+Dfe(2*n);
- b) a# i/ X$ H+ E  j5 n+ |    end
    f
  P3 s; a, O  c" j3 D) }8 @    e
1 L7 w& t" B7 {% aend


8 G1 Z: O  [0 F: P9 o运行会出现
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

% V, j  w" l1 I2 H0 N8 [6 o! W  j我这个程序是按照陈珩那本稳态分析的例子编的，想试试。完全没考虑PV节点
$ e. }0 {. y; b% C6 o这个程序要怎么改进啊？

wangmuyidong

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