matlab实现潮流计算程序.txt

这是一段matlab程序

Matlab实现潮流计算程序你可以通过这个链接引用该篇文章:http://lntianxia.bokee.com/viewdiary.16330339.html

程序代码如下:

111111.

 %读入数据
clc
clear
filename='123.txt';
a=textread(filename)

n=a(1,1);
pinghengjd=a(1,2);
phjddianya=a(1,3);
jingdu=a(1,4);

b=zeros(1,9);
j1=0;
[m1,n1]=size(a);
for i1=1:m1
    if a(i1,1)==0
        j1=j1+1;        
        b(j1)=i1;
    end
end
b;

%矩阵分块
a1=a(b(1)+1:b(2)-b(1)+1,1:n1);
a2=a(b(2)+1:b(3)-1,1:n1);
a3=a(b(3)+1:b(4)-1,1:n1);
a4=a(b(4)+1:b(5)-1,1:n1);
a5=a(b(5)+1:b(6)-1,1:n1);

%设置初值
vcz=1;
dcz=0;
kmax=20;
k1=0;

%求节点导纳矩阵
a11=zeros(4,6);
for i0=1:3
    for j0=1:6
        a11(i0,j0)=a1(i0,j0);
        a11(4,j0)=a2(1,j0);
    end
end
a11;

linei=a11(1:4,2);
linej=a11(1:4,3);
liner=a11(1:4,4);
linex=a11(1:4,5);
lineb=a11(1:4,6);

branchi=0;
branchj=0;
branchb=0;

G=zeros(4,4);
B=zeros(4,4);

for k=1:4
    i2=linei(k,1);
    j2=linej(k,1);
    r=liner(k,1);
    x=linex(k,1);
    b=0;
    GIJ=r/(r*r+x*x);
    BIJ=-x/(r*r+x*x);
    if k>=4 & lineb(k)~=0
        k0=lineb(k);
        G(i2,j2)=-GIJ/k0;
        G(j2,i2)=G(i2,j2);
        B(i2,j2)=-BIJ/k0;
        B(j2,i2)=B(i2,j2);
        G(i2,i2)=G(i2,i2)+GIJ/k0/k0;
        B(i2,i2)=B(i2,i2)+BIJ/k0/k0;
    else
        G(j2,i2)=-GIJ;
        G(i2,j2)=G(j2,i2);
        B(j2,i2)=-BIJ;
        B(i2,j2)=B(j2,i2);
        G(i2,i2)=G(i2,i2)+GIJ;
        b=lineb(k);
        B(i2,i2)=B(i2,i2)+BIJ+b;
    end
    G(j2,j2)=G(j2,j2)+GIJ;
    B(j2,j2)=B(j2,j2)+BIJ+b;
end
G;
B;

B=B.*i;
Yf=G+B
Y=abs(Yf);
alf=angle(Yf);

%赋Jacobian矩阵参数
P=zeros(n,1);
Q=zeros(n,1);

Pd=zeros(1,n);
Qd=zeros(1,n);
dP=zeros(1,n);
dQ=zeros(1,n);

PG=a4(:,3);
PD=a4(:,5);
QG=a4(:,4);
QD=a4(:,6);
i8=a4(:,2);

for j8=1:length(i8)
    P(i8(j8))=PG(i8(j8))-PD(i8(j8));
    Q(i8(j8))=QG(i8(j8))-QD(i8(j8));
end

delt=zeros(n,1);
V=ones(n,1);
V(3)=1.10;
V(4)=1.05;

ddelt=zeros(n,1);
dV=zeros(n,1);
A=zeros(2*n,2*n);
B=zeros(2*n,1);

Jacobian=Jaco(V,delt,n,Y,alf)

%求取矩阵功率
for j5=1:kmax
    disp(['第' int2str(j5) '次计算结果'])   
    if k>=kmax
        break
    end    
    for i10=1:4
        Pd(i10)=0;
        Qd(i10)=0;
        for j10=1:n
            Pd(i10)=Pd(i10)+V(i10)*Y(i10,j10)*V(j10)*cos(delt(i10)-delt(j10)-alf(i10,j10));
            Qd(i10)=Qd(i10)+V(i10)*Y(i10,j10)*V(j10)*sin(delt(i10)-delt(j10)-alf(i10,j10));                
        end
    end
    for i4=1:3
        dP(i4)=P(i4)-Pd(i4);            
    end
    for j4=1:2
        dQ(j4)=Q(j4)-Qd(j4);
    end
    A=Jaco(V,delt,n,Y,alf)
    for i14=1:n
        B(i14*2-1)=-dP(i14);
        B(i14*2)=-dQ(i14);
    end
    if max(abs(B))>jingdu
        X=A\B;
        for i16=1:n
            ddelt(i16)=X(2*i16-1);
            dV(i16)=X(2*i16)*V(i16);
        end
         V=V+dV
         delt=delt+ddelt 
    else
         break
    end
    disp('----------------')
end   
    
    
    
%流氓算法        
%        for ii=1:2  
%            V(ii)=V(ii)+dV(ii);        
%        end   
%        V   

 

222222.

function A=Jaco(V,delt,n,Y,alf)
%计算Jacobian矩阵
for i7=1:n
    Hd1(i7)=0;
    Jd1(i7)=0;
    for j7=1:n
        Hd1(i7)=Hd1(i7)+V(i7)*Y(i7,j7)*V(j7)*sin(delt(i7)-delt(j7)-alf(i7,j7));
        Jd1(i7)=Jd1(i7)+V(i7)*Y(i7,j7)*V(j7)*cos(delt(i7)-delt(j7)-alf(i7,j7));
    end
end

for i6=1:n
    for j6=1:n
        if i6~=j6
            H(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*sin(delt(i6)-delt(j6)-alf(i6,j6));
            N(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6));
            J(i6,j6)=-N(i6,j6);
            L(i6,j6)=H(i6,j6);
        else
            H(i6,i6)=Hd1(i6)-V(i6)*Y(i6,i6)*V(i6)*sin(delt(i6)-delt(j6)-alf(i6,j6));
            J(i6,j6)=-Jd1(i6)+V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6));
            N(i6,j6)=-Jd1(i6)-V(i6)*Y(i6,i6)*V(i6)*cos(alf(i6,i6));
            L(i6,i6)=-Hd1(i6)+V(i6)*Y(i6,i6)*V(i6)*sin(alf(i6,i6));
        end
    end
end

%修正Jacobian矩阵
for j9=3
    for i9=1:n
        N(i9,j9)=0;
        L(i9,j9)=0;
        J(j9,i9)=0;
        L(j9,i9)=0;
    end
end
L(j9,j9)=1;

for j9=4
    for i9=1:n
        H(i9,j9)=0;
        N(i9,j9)=0;
        J(i9,j9)=0;
        L(i9,j9)=0;
        H(j9,i9)=0;
        N(j9,i9)=0;
        J(j9,i9)=0;
        L(j9,i9)=0;
    end
end
H(j9,j9)=1;
L(j9,j9)=1;
%Jaco=[H N;J L];
%Jaco=zeros(2*n,2*n);
for i11=1:n
    for j11=1:n
        Jaco(2*i11-1,2*j11-1)=H(i11,j11);
        Jaco(2*i11-1,2*j11)=N(i11,j11);
        Jaco(2*i11,2*j11-1)=J(i11,j11);
        Jaco(2*i11,2*j11)=L(i11,j11);
    end
end
A=Jaco;

33333.

数据:

4      4      1.05    0.00001
0
1      1      2       0.1      0.40     0.01528
2      1      4       0.12     0.50     0.01920 
3      2      4       0.08     0.40     0.01413
0
1      1      3       0        0.3      0.90909091
0
0
1      1      0       0        0.30     0.18
2      2      0       0        0.55     0.13
3      3      0.5     0        0        0
0
1      3      1.10    0        0
0