clc
clear all
fid=fopen('rts.txt','r');
N=fscanf(fid,'%f',[1,1]);
NL=fscanf(fid,'%f',[1,1]);
AA=fscanf(fid,'%f',[9,N]);SS=AA';
BB=fscanf(fid,'%f',[10,NL]);ZL=BB';
fclose(fid);
MM=zeros(N,1);NN=zeros(N,1);
cc=0;
N2=N;
NL2=NL;
E=zeros(1,10);
fz=[];
xj=[];
for i=1:N
    if(SS(i,3)~=0||SS(i,4)~=0)
        cc=cc+1;
        N2=N2+1;
        NL2=NL2+1; 
        E(1,cc)=SS(i,7)+j*SS(i,9)*(SS(i,3)-SS(i,4)*j)/SS(i,7);
        fz(1,cc)=abs(E(1,cc));
        xj(1,cc)=angle(E(1,cc))*(180/pi);
        SS(N2,1)=N2;
        SS(N2,2)=1;
        SS(N2,3)=0;
        SS(N2,4)=0;
        SS(N2,5)=0;
        SS(N2,6)=0;
        SS(N2,7)=fz(1,cc);
        SS(N2,8)=xj(1,cc); % 發電機內電勢.
        SS(N2,9)=SS(i,9);
        ZL(NL2,1)=i;
        ZL(NL2,2)=N2;
        ZL(NL2,3)=0;
        ZL(NL2,4)=SS(N2,9); % 暫態電抗
        ZL(NL2,5)=0; % 5-10 自行輸入.
        ZL(NL2,6)=0;
        ZL(NL2,7)=0;
        ZL(NL2,8)=0;
        ZL(NL2,9)=0;
        ZL(NL2,10)=0;
    end    
end
N=N2;
NL=NL2; 
SS
ZL
Y=zeros(N);
for t=1:NL    %線路導納
    I=ZL(t,1);
    J=ZL(t,2);
    Y(I,I)=Y(I,I)+1/(ZL(t,3)+(ZL(t,4))*j)+j*0.5*ZL(t,5);
    Y(J,J)=Y(J,J)+1/(ZL(t,3)+(ZL(t,4))*j)+j*0.5*ZL(t,5);
    Y(I,J)=-1/(ZL(t,3)+ZL(t,4)*j);
    Y(J,I)=-1/(ZL(t,3)+ZL(t,4)*j);
end
YGG=Y(25:34,25:34);
YGL=Y(1:24,25:34);
YLG=Y(25:34,1:24);
YLL=Y(1:24,1:24);
Y5=YGG-YLG*inv(YLL)*YGL;
Y5
G=real(Y5);
B=imag(Y5);
Kii=zeros(1,10);
for i=1:10
    for j=1:10
       Kii(1,i)=Kii(1,i)+((-G(i,j)*sin(xj(1,i)-xj(1,j))+B(i,j)*cos(xj(1,i)-xj(1,j)))*fz(1,j)-B(i,i)*fz(1,i))*fz(1,i);
    end
end
Kij=zeros(10);
for i=1:10
    for j=1:10
        if i~=j
            Kij(i,j)=fz(1,i)*fz(1,j)*(G(i,j)*sin(xj(1,i)-xj(1,j))-B(i,j)*cos(xj(1,i)-xj(1,j)));
        else
            Kij(i,j)=Kii(1,i);
        end
    end
end
Sb=1.00;
D=[2,2,2,2,2,2,2,2,2,2];
Wk=[23.50,117.00,183.00,254.00,107.90,168.00,176.00,524.00,233.00,105.00];
Tj=2*Wk/Sb;
C=zeros(20); 
C1=zeros(10);
C2=eye(10);
C3=zeros(10);
C4=eye(10);
for i=1:10
    for j=1:10
        C3(i,j)=-Kij(i,j)/Tj(i);
    end
end
for n=1:10
    C4(n,n)=-D(n)/Tj(n);
end
C=[C1 C2;C3 C4];
A=eig(C);
A
for i=1:10
    plot(A(i),'b*');
    hold on;
end
CC=[-0.025,0;0.002,0];plot(CC(:,1),CC(:,2),'r');
CC=[0,-1;0,1];plot(CC(:,1),CC(:,2),'k');
box off;
axis([-0.025,0.002,-1,1]);
xlabel('實部');
ylabel('虛部');
title('電力系統靜態穩定特征值')
            
break;
bar(A);