azhen.m

对多机电力系统进行小干扰稳定分析 可以计算特征值

function Azhen
clear;
clc;
b=[4 5 0.01+0.085i 0.176i 1;
    4 6 0.017+0.092i 0.158i 1;
    5 7 0.032+0.161i 0.306i 1;
    6 9 0.039+0.17i 0.358i 1;
    7 8 0.0085+0.072i 0.149i 1;
    8 9 0.0119+0.1008i 0.209i 1;
    1 4 0.0576i 0 1;
    2 7 0.0625i 0 1;
    3 9 0.0586i 0 1];
%——b中数据依次为——支路始端号、支路末端号、支路阻抗、线路电容、支路变比
nl=size(b,1);
n=9;%系统节点数
n1=3;%发电机数
%disp(nl)
%disp(n)
Y=zeros(n);
for i=1:nl%导纳矩阵生成
    p=b(i,1);q=b(i,2);
    Y(p,q)=Y(p,q)-1./(b(i,3)*b(i,5));
    Y(q,p)=Y(p,q);
    Y(q,q)=Y(q,q)+1./b(i,3)+b(i,4)./2;
    Y(p,p)=Y(p,p)+1./(b(i,3)*b(i,5)^2)+b(i,4)./2;
end
s1=[1.25+0.5i;0.9+0.3i;1+0.35i];%负荷功率
u1=[0.99563;1.01265;1.01588];%负荷节点电压幅值
Y(5,5)=Y(5,5)+conj(s1(1))./(u1(1)^2);
Y(6,6)=Y(6,6)+conj(s1(2))./(u1(2)^2);
Y(8,8)=Y(8,8)+conj(s1(3))./(u1(3)^2);

disp('并入负荷节点后的Y阵：')
disp(Y)
G=real(Y);B=imag(Y);
%disp(G)
for i=1:n
    for j=1:n
        Y1(2*i-1,2*j-1)=G(i,j);
        Y1(2*i-1,2*j)=-B(i,j);
        Y1(2*i,2*j-1)=B(i,j);
        Y1(2*i,2*j)=G(i,j);
    end
end
disp('将负荷节点并入Y阵后增阶的Y阵：')
disp(Y1)
    for i=1:2*n1%发电机节点
        for j=1:2*n1
            Ynn(i,j)=Y1(i,j);
        end
    end
    for i=1:2*n-2*n1%剩余节点1
        for j=1:2*n-2*n1
            Yrr(i,j)=Y1(i+2*n1,j+2*n1);
        end
    end
    for i=1:2*n1%剩余节点2
        for j=1:2*n-2*n1
            Ynr(i,j)=Y1(i,j+2*n1);
        end
    end
    for i=1:2*n-2*n1%剩余节点3
        for j=1:2*n1
            Yrn(i,j)=Y1(i+2*n1,j);
        end
    end
    ooooo=size(Yrr)
    disp('发电机节点：')
    disp(Ynn)
    disp('剩余节点1：')
    disp(Yrr)
    disp('剩余节点2：')
    disp(Ynr)
    disp('剩余节点3：')
    disp(Yrn)
Yx=Ynn-Ynr*inv(Yrr)*Yrn;
disp('消去所有联络节点后的Yx阵：')
disp(Yx)
d=[3.58572*pi./180;61.09844*pi./180;54.13662*pi./180];
%d=[2.27165*pi./180;61.09844*pi./180;54.13662*pi./180];
%d=[2.27165*pi./180;19.73159*pi./180;13.16641*pi./180];
xd1=[0.0608;0.1198;0.1813];
xq=[0.0969;0.8645;1.2578];
for i=1:n1
    gf1(i)=-(xd1(i)-xq(i))*sin(2*d(i))./(2*xd1(i)*xq(i));
    bf1(i)=(xd1(i)*(cos(d(i))^2)+xq(i)*(sin(d(i))^2))./(xd1(i)*xq(i));%此时已将bf1变为-bf1
    bf2(i)=-(xd1(i)*(sin(d(i))^2)+xq(i)*(cos(d(i))^2))./(xd1(i)*xq(i));
    gf2(i)=(xd1(i)-xq(i))*sin(2*d(i))./(2*xd1(i)*xq(i));
end
% disp(gf1)
% disp(bf1)
% disp(bf2)
% disp(gf2)
for i=1:n1%将发电机参数并入导纳阵
    Yx(2*i-1,2*i-1)=Yx(2*i-1,2*i-1)+gf1(i);
    Yx(2*i-1,2*i)=Yx(2*i-1,2*i)+bf1(i);
    Yx(2*i,2*i-1)=Yx(2*i,2*i-1)+bf2(i);
    Yx(2*i,2*i)=Yx(2*i,2*i)+gf2(i);
end
disp('将发电机参数并入导纳阵后的Yx阵：')
disp(Yx)
    Yn=inv(Yx);
    disp('导纳阵的逆阵Yn：')
    disp(Yn)
    for i=1:n1%所求为gf1等对功角的倒数
        gf11(i)=-(xd1(i)-xq(i))*cos(2*d(i))./(xd1(i)*xq(i));
    bf11(i)=(xd1(i)-xq(i))*sin(2*d(i))./(xd1(i)*xq(i));
    bf21(i)=-(xd1(i)-xq(i))*sin(2*d(i))./(xd1(i)*xq(i));
    gf21(i)=(xd1(i)-xq(i))*cos(2*d(i))./(xd1(i)*xq(i));%分母多了个2
    end
%     disp(gf11)
% disp(bf11)
% disp(bf21)
% disp(gf21)
for i=1:n1%由dq变换到xy的坐标变换矩阵
    T(2*i-1,2*i-1)=sin(d(i));
    T(2*i-1,2*i)=cos(d(i));
    T(2*i,2*i-1)=-cos(d(i));
    T(2*i,2*i)=sin(d(i));
end
% disp('变换矩阵：')
% disp(T)
T1=inv(T);

ux=zeros(n1,1);
uy=zeros(n1,1);
uu=zeros(n1,1);
Ii=zeros(n1,1);
Uxy=zeros(2*n1,1);
Ixy=zeros(2*n1,1);
Udq=zeros(2*n1,1);
Idq=zeros(2*n1,1);
d1=[0;9.28*pi./180;4.6648*pi./180];%发电机端电压相角
u=[1.04;1.025;1.025];%发电机端电压幅值
for i=1:n1
    ux(i)=u(i)*cos(d1(i));
    uy(i)=u(i)*sin(d1(i));
end
% disp('ux')
% disp(ux)
% disp('uy')
% disp(uy)
s=[0.71641+0.27046i;1.63+0.06654i;0.85-0.1086i];%发电机功率
uu=complex(ux,uy);
% disp('uu')
% disp(uu)
for i=1:n1
    Ii(i)=conj(s(i))./conj(uu(i));
end
% disp('Ii')
% disp(Ii)
ix=real(Ii);
iy=imag(Ii);
% disp(ix)
% disp(iy)
for i=1:n1
Uxy(2*i-1)=ux(i);
Uxy(2*i)=uy(i);
Ixy(2*i-1)=ix(i);
Ixy(2*i)=iy(i);
end
for i=1:n1
    Udq=T1*Uxy;
    Idq=T1*Ixy;
end
% disp('Udq')
% disp(Udq)
% disp('Idq')
% disp(Idq)
for i=1:n1
    uq(i)=Udq(2*i);
    id(i)=Idq(2*i-1);
end
% disp('uq')
% disp(uq)
% disp('id')
% disp(id)
for i=1:n1
    eq1(i)=uq(i)+xd1(i)*id(i);
end
% disp('eq1')
% disp(eq1)

for i=1:n1
    aq(i)=gf1(i)*cos(d(i))+bf1(i)*sin(d(i));%bf1前系数应该为正，因为前面所求的bf1等于-bf1   
    bq(i)=bf2(i)*cos(d(i))+gf2(i)*sin(d(i));
    ad(i)=eq1(i)*[(gf11(i)+bf1(i))*cos(d(i))-(bf11(i)+gf1(i))*sin(d(i))]-(gf11(i)*ux(i)-bf11(i)*uy(i));
    bd(i)=eq1(i)*[(bf21(i)+gf2(i))*cos(d(i))+(gf21(i)-bf2(i))*sin(d(i))]-(bf21(i)*ux(i)+gf21(i)*uy(i));
end
% disp('参数aq:')
% disp(aq)
% disp('参数bq:')
% disp(bq)
disp('参数ad:')
disp(ad)
disp('参数bd:')
disp(bd)
for i=1:n1
    for j=1:n1
        cq(i,j)=Yn(2*i-1,2*j-1)*aq(j)+Yn(2*i-1,2*j)*bq(j);
        dq(i,j)=Yn(2*i,2*j-1)*aq(j)+Yn(2*i,2*j)*bq(j);
        cd(i,j)=Yn(2*i-1,2*j-1)*ad(j)+Yn(2*i-1,2*j)*bd(j);
        dd(i,j)=Yn(2*i,2*j-1)*ad(j)+Yn(2*i,2*j)*bd(j);
    end
end
% disp('参数cq:')
% disp(cq)
% disp('参数dq:')
% disp(dq)
disp('参数cd:')
disp(cd)
disp('参数dd:')
disp(dd)
for i=1:n1
   for j=1:n1
       if i~=j
        eq(i,j)=-(gf1(i)*cq(i,j)+bf1(i)*dq(i,j));
        fq(i,j)=-(bf2(i)*cq(i,j)+gf2(i)*dq(i,j));
        ed(i,j)=-(gf1(i)*cd(i,j)+bf1(i)*dd(i,j));
        fd(i,j)=-(bf2(i)*cd(i,j)+gf2(i)*dd(i,j));
       else
    eq(i,j)=-(gf1(i)*cq(i,i)+bf1(i)*dq(i,i))+aq(i);
    fq(i,j)=-(bf2(i)*cq(i,i)+gf2(i)*dq(i,i))+bq(i);
    ed(i,j)=-(gf1(i)*cd(i,i)+bf1(i)*dd(i,i))+ad(i);
    fd(i,j)=-(bf2(i)*cd(i,i)+gf2(i)*dd(i,i))+bd(i);
       end
    end
end
% disp('参数eq:')
% disp(eq)
% disp('参数fq:')
% disp(fq)
disp('参数ed:')
disp(ed)
disp('参数fd:')
disp(fd)
for i=1:n1%求系数k1所需的diag(Uxi)等矩阵
    ux1(i,i)=ux(i);
    uy1(i,i)=uy(i);
    ix1(i,i)=ix(i);
    iy1(i,i)=iy(i);
end
% disp('求系数k1所需的diag(Uxi)等矩阵')
% disp(ux1)
% disp(uy1)
% disp(ix1)
% disp(iy1)
k1=ux1*ed+uy1*fd+ix1*cd+iy1*dd;%每一行元素相加所得的值为零
k2=ux1*eq+uy1*fq+ix1*cq+iy1*dq;
disp('系数k1:')
disp(k1)
for i=1:n1
    a=0;
    for j=1:n1
        a=a+k1(i,j);
    end
    disp('k1每行相加结果：')
disp(a)
end
disp('系数k2:')
disp(k2)
xd=[0.146;0.8958;1.3125];
for i=1:n1
    dk(i,i)=xd(i)-xd1(i);%此式为系数k4公式中的第一项
    cs(i,i)=cos(d(i))*ix(i)+sin(d(i))*iy(i);%此式为系数k4公式中的第二项
    ds(i,i)=sin(d(i));%此式为系数k4公式中的第三项
    dc(i,i)=cos(d(i));%此式为系数k4公式中的第四项
end
%disp(dk)
k4=dk*(cs+ds*ed-dc*fd);%每一行元素相加所得的值为零
% disp('系数k4:')
% disp(k4)
% for i=1:n1
%     aa=0;
%     for j=1:n1
%         aa=aa+k4(i,j);
%     end
%     disp('k4每行相加结果：')
% disp(aa)
% end
k3=eye(n1)+dk*(ds*eq-dc*fq);
% k3=inv(k3);
% disp('系数k3:')
% disp(k3)
for i=1:n1
    ut(i)=(ux(i)^2+uy(i)^2)^0.5;
end
disp('发电机端电压Ut：')
disp(ut)
for i=1:n1
    du(i,i)=ux(i)./ut(i);%k5表达式中第一项
    du1(i,i)=uy(i)./ut(i);%k5表达式中第二项
end
% disp(du)
k5=du*cd+du1*dd;
k6=du*cq+du1*dq;
% disp('系数k5:')%每一行元素相加所得的值为零
% disp(k5)
% for i=1:n1
%     aaa=0;
%     for j=1:n1
%         aaa=aaa+k5(i,j);
%     end
%     disp('k5每行相加结果：')
% disp(aaa)
% end
% disp('系数k6:')
% disp(k6)
tj=[47.28;12.8;6.02];
td1=[8.96;6;5.89];
% tq1=[0;0.535;0.6];
D=[0;0;0];
%D=[1;1;1];
% ka=[50;50;50];
ka=[500;500;500];
ta=[0.03;0.03;0.03];

te=[0.5;0.5;0.5];
kf=[0.04;0.04;0.04];
tf=[0.715;0.715;0.715];

A13=zeros(n1);D11=zeros(n1);
for i=1:n1
    M(i,i)=tj(i);
    D1(i,i)=D(i);%式11-52中的D
    Td01(i,i)=td1(i);
    A11(i,i)=-1./ta(i);
    A13(i,i)=-ka(i)./ta(i);
    A21(i,i)=1./te(i);
    A31(i,i)=kf(i)./(tf(i)*te(i));
    A33(i,i)=-1./tf(i);
end
    D11=A13*k5;
    D12=zeros(n1);
    D13=A13*k6;
    A32=-A31;
    A22=-A21;
%     disp('M')
%     disp(M)
disp('A阵中参数：')
% disp(D1)
% disp(Td01)
% disp(A11)
% disp(A13)
% disp(A21)
% disp(A22)
% disp('A31:')
% disp(A31)
% disp('A32:')
% disp(A32)
% disp(A33)
% disp(D11)
% disp(D12)
% disp(D13)
Mk1=-inv(M)*k1;
% disp('Mk1')
% disp(Mk1)
MD=-inv(M)*D1;
Mk2=-inv(M)*k2;
Td01k4=-inv(Td01)*k4;%应该是-Td01乘以k4和k3
Td01k3=-inv(Td01)*k3;
Td0=inv(Td01);
f0=50;
I=zeros(n1);
I1=[2*pi*f0,2*pi*f0,2*pi*f0];
I=diag(I1);
% disp(I)

A=zeros(6*n1);%A阵的生成
for i=1:n1
    for j=n1+1:1:n1+n1
        A(i,j)=I(i,j-n1);
    end
end
%disp(A)
for i=n1+1:1:n1+n1
    for j=1:n1
        A(i,j)=Mk1(i-n1,j);
    end
end
%disp(A)
for i=n1+1:1:n1+n1
    for j=n1+1:1:n1+n1
        A(i,j)=MD(i-n1,j-n1);
    end
end
for i=n1+1:1:n1+n1
    for j=2*n1+1:1:2*n1+n1
        A(i,j)=Mk2(i-n1,j-2*n1);
    end
end

for i=2*n1+1:1:2*n1+n1
    for j=1:n1
        A(i,j)=Td01k4(i-2*n1,j);
    end
end

for i=2*n1+1:1:2*n1+n1
    for j=2*n1+1:1:2*n1+n1
    A(i,j)=Td01k3(i-2*n1,j-2*n1);
    end
end

for i=3*n1+1:1:3*n1+n1
    for j=1:n1
        A(i,j)=D11(i-3*n1,j);
    end
end

for i=3*n1+1:1:3*n1+n1
    for j=n1+1:1:n1+n1
        A(i,j)=D12(i-3*n1,j-n1);
    end
end

for i=3*n1+1:1:3*n1+n1
    for j=2*n1+1:1:2*n1+n1
        A(i,j)=D13(i-3*n1,j-2*n1);
    end
end

for i=2*n1+1:1:2*n1+n1
    for j=4*n1+1:1:4*n1+n1
        A(i,j)=Td0(i-2*n1,j-4*n1);
    end
end

for i=3*n1+1:1:3*n1+n1
    for j=3*n1+1:1:3*n1+n1
        A(i,j)=A11(i-3*n1,j-3*n1);
    end
end

for i=3*n1+1:1:3*n1+n1
    for j=5*n1+1:1:5*n1+n1
        A(i,j)=A13(i-3*n1,j-5*n1);
    end
end

for i=4*n1+1:1:4*n1+n1
    for j=3*n1+1:1:3*n1+n1
        A(i,j)=A21(i-4*n1,j-3*n1);
    end
end

for i=4*n1+1:1:4*n1+n1
    for j=4*n1+1:1:4*n1+n1
        A(i,j)=A22(i-4*n1,j-4*n1);
    end
end

for i=5*n1+1:1:5*n1+n1
    for j=3*n1+1:1:3*n1+n1
        A(i,j)=A31(i-5*n1,j-3*n1);
    end
end

for i=5*n1+1:1:5*n1+n1
    for j=4*n1+1:1:4*n1+n1
        A(i,j)=A32(i-5*n1,j-4*n1);
    end
end

for i=5*n1+1:1:5*n1+n1
    for j=5*n1+1:1:5*n1+n1
        A(i,j)=A33(i-5*n1,j-5*n1);
    end
end
% disp('整个系统的状态矩阵A：')
% disp(A)
% disp('6阶全部特征值：')
eig(A)
for i=1:2*n1
    for j=1:2*n1
        A1(i,j)=A(i,j);
    end
end
% disp('经典模型的状态矩阵：')
% disp(A1)
% disp('部分特征值：')
eig(A1)