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

zhouyou529

clear;
! d5 O8 m$ q& p$ o& k% u4 U: I. xYb=[6.25-18.75i  -5+15i        -1.25+3.75i    0             0
( _  y( y7 v+ D7 n& S    -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
3 R  c$ E; G$ p6 n2 ^' |    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];
Pknow=[0;0.2;-0.45;-0.4;-0.6];Qknow=[0;0.2;-0.15;-0.05;-0.1];e1=1.06;f1=0;%已知量
) v- }" G# ]3 N8 z0 c4 [e=[1.06;1;1;1;1];f=[0;0;0;0;0];
" N4 F/ b/ z" O' L' M" Ik=0;%跌代次数
. T) M4 \7 D+ Y3 \4 VPQX=[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;
, h& t5 Q! _5 c% G+ w1 m/ Q& H    for n=1:4
        for m=1:5
- J4 Q4 Z4 n. u& I) X& C6 W4 O            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)));
. n" B7 \7 i  O        end
        PQX(2*n-1)=Pknow(n+1)-P(n);
        PQX(2*n)=Qknow(n+1)-Q(n);
        I(n)=(P(n)-Q(n)*i)/(e(n+1)-f(n+1)*i);
    end
$ Y( }, T; c2 X7 ?" e. c    for n=1:4
( T6 I. F0 V# L; y( s        for m=1:4
, c% ?# |3 Z4 f4 v8 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
* K3 c; K9 X3 I% L2 M            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
0 y' C0 P, l4 r            YAKEBI(2*n,2*m)=-imag(Yb(n+1,m+1))*e(n+1)+real(Yb(n+1,m+1))*f(n+1);%Lnm
- t% O/ |" _: T        end
* A$ \7 E; ~( ~" T# j0 Y7 m/ m        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
& {6 B) @; r# Y( Q! c! o        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
$ |: p, A) Q# b! m+ v1 E, L' W    Dfe=inv(YAKEBI)*PQX; %修正fe
; m% ]' g8 c0 j2 U8 l( j+ R( ?    for n=1:4
2 o8 n$ c) ^8 @% m7 i( }5 l; C3 O$ q        f(n+1)=f(n+1)+Dfe(2*n-1);
        e(n+1)=e(n+1)+Dfe(2*n);
    end
    f
7 z7 [' p) s/ {6 N9 Z) y    e
& ?2 e. ~* f; L  ~+ q- i7 Uend

5 {9 u4 R& k  d7 q& d: T3 D5 w" D
运行会出现
3 r) d5 U' s% W+ [, n' |  F1 v7 cWarning: Matrix is singular to working precision.
, k2 l4 k) U: X" F( e# u; H* k. a(Type "warning off MATLAB:singularMatrix" to suppress this warning.)
  D, j& D% ]% l& L3 M$ b> In f:\MATLAB6p5\work\yy.m at line 35
1 r3 ]+ {7 b' Z3 g/ G+ |
我这个程序是按照陈珩那本稳态分析的例子编的，想试试。完全没考虑PV节点
这个程序要怎么改进啊？

wangmuyidong

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