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

zhouyou529

clear;
% z* g6 M( @+ ^4 D: EYb=[6.25-18.75i  -5+15i        -1.25+3.75i    0             0
7 H0 P8 R+ @# M4 e0 Z    -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
    0            -1.667+5i     -10+30i       12.917-38.75i  -1.25+3.75i
& B9 h6 [1 X; _+ s. h9 T5 a    0            -2.5+7.5i     0             -1.25+3.75i    3.75-11.25i];
2 _. r6 D: B) e; 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;%已知量
* H% |/ y5 F- f3 D2 V5 ze=[1.06;1;1;1;1];f=[0;0;0;0;0];
k=0;%跌代次数
( `. q6 M% u9 c& r6 J0 EPQX=[1;1;1;1;1;1;1;1];%修正量   
; |/ R2 l: E9 OP=zeros(4,1);Q=zeros(4,1);I=zeros(4,1);YAKEBI=zeros(8,8);
while (max(PQX)>=0.1|min(PQX)<=-0.1)
2 O" {1 Z& H: P' u' {/ M  c    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)));
1 U" [  t5 j( M" {  a- X            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)));
: q- B+ D9 }* g9 W% P        end
# r. ~7 e  O3 p) a8 P/ L( W/ x4 w        PQX(2*n-1)=Pknow(n+1)-P(n);
: w5 p. A' H" Y0 P" v3 [0 }        PQX(2*n)=Qknow(n+1)-Q(n);
" F  r6 A) ]7 M2 u        I(n)=(P(n)-Q(n)*i)/(e(n+1)-f(n+1)*i);
    end
    for n=1:4
        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
5 {. L8 H+ `  d# 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
" l% ~0 K$ [  z& q0 v            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
6 B  K- k2 H8 M( Y' 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
: n# o+ k' J# [: l& S/ F4 k+ u: f2 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
& J7 F1 R1 ?/ D2 K7 U        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
; _: ]) F0 C  d- F    end% build the yakebi matrix
; l; F6 I) S/ E' p- k. V$ C3 E    Dfe=inv(YAKEBI)*PQX; %修正fe
# S7 R+ _( J1 |. l! {; M/ |& y    for n=1:4
; b( u" z% I& M: ^        f(n+1)=f(n+1)+Dfe(2*n-1);
( N! H8 s. J( ?        e(n+1)=e(n+1)+Dfe(2*n);
    end
7 I) n  V: [3 b$ F3 {7 G0 f  u    f
0 r  M; W* t) t  t( A4 u8 R    e
end
% B: v3 H/ S" H4 m' L7 A

! B2 n$ H8 Q( e& @9 H运行会出现
Warning: Matrix is singular to working precision.
# Z5 \! [$ Q7 ~2 L(Type "warning off MATLAB:singularMatrix" to suppress this warning.)
> In f:\MATLAB6p5\work\yy.m at line 35
: I2 Z3 h7 n& `1 z. V
我这个程序是按照陈珩那本稳态分析的例子编的，想试试。完全没考虑PV节点
+ b, ]8 l, r3 K1 z这个程序要怎么改进啊？

wangmuyidong

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