ncewscc2.m

可以计算电力系统概率潮流

clc
clear

clear global
% global variables
glo_var
% Edgeworth's Expansion Theory
clc
nbus=179;
nline=222;
t=cputime;

bus=csvread('busdat.csv');
% bus(:,1) column 1 is the system index
% bus(:,8) column 8 is the program index

line1=csvread('linedat.csv');
% modline is to simplify the line reactance of multiple lines to one 
line=modline(line1);


% find Swing Bus Position
% nswing is the swing bus number
nswing=findSwing(bus);

% Y is the admittance matrix
Y=FormY(bus, line);
% redY is the the reduced admittance matrix without swing bus column and row
Yred=redY(Y, nswing);
Gy=real(Y);
By=imag(Y);

% bus zhong de bus(:,3) wei jiaodu,yao zhuan huan wei hu du zhi.
bus(:,3)=bus(:,3)*pi/180;

% gen ju Qinghua P142 
% qiu zhu ru gonglv de yakebi juzhen
H1(179,179)=0;
N1(179,179)=0;
J1(179,179)=0;
L1(179,179)=0;
for i=1:179
    for j=1:179
        if i==j
            H1(i,i)=bus(i,2)^2*By(i,i)+bus(i,7)-bus(i,5);
            N1(i,i)=-bus(i,2)^2*Gy(i,i)-bus(i,6)+bus(i,4);
            J1(i,i)=bus(i,2)^2*Gy(i,i)-bus(i,6)+bus(i,4);
            L1(i,i)=bus(i,2)^2*By(i,i)-bus(i,7)+bus(i,5);
        else
            H1(i,j)=-bus(i,2)*bus(j,2)*(Gy(i,j)*sin(bus(i,3)-bus(j,3))-By(i,j)*cos(bus(i,3)-bus(j,3)));
            N1(i,j)=-bus(i,2)*bus(j,2)*(Gy(i,j)*cos(bus(i,3)-bus(j,3))+By(i,j)*sin(bus(i,3)-bus(j,3)));
            J1(i,j)=bus(i,2)*bus(j,2)*(Gy(i,j)*cos(bus(i,3)-bus(j,3))+By(i,j)*sin(bus(i,3)-bus(j,3)));
            L1(i,j)=-bus(i,2)*bus(j,2)*(Gy(i,j)*sin(bus(i,3)-bus(j,3))-By(i,j)*cos(bus(i,3)-bus(j,3)));
        end
    end
end
H1=-H1;
N1=-N1;
J1=-J1;
L1=-L1;
H1=redJ0(H1, nswing);
N1=redJ0(N1, nswing);
J1=redJ0(J1, nswing);
L1=redJ0(L1, nswing);
% Yakebi da Juzhen
Ykb(1:178,1:178)=H1;
Ykb(179:356,1:178)=J1;
Ykb(1:178,179:356)=N1;
Ykb(179:356,179:356)=L1;
% Ykb de ni juzhen
S0=inv(Ykb);

% zhilu guan xi juzhen

H2(178,178)=0;
N2(178,178)=0;
J2(178,178)=0;
L2(178,178)=0;
for i=1:nline
    frompos=line(i,1);
    topos=line(i,2);
    startpos=FindPos(bus, frompos);
    endpos=FindPos(bus, topos);
    if (startpos>nswing)
        H2(i,startpos-1)=bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3))-By(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3)));
        N2(i,startpos-1)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3))+By(startpos,endpos)*sin(bus(startpos,3)-bus(j,3)))+2*bus(startpos,2)^2*By(startpos,endpos);
        J2(i,startpos-1)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3))+By(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3)));
        L2(i,startpos-1)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3))-By(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3)))-2*bus(startpos,2)^2*line(i,7)+bus(startpos,2)^2*By(startpos,endpos);
        
    elseif (startpos<nswing)
        H2(i,startpos)=bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3))-By(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3)));
        N2(i,startpos)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3))+By(startpos,endpos)*sin(bus(startpos,3)-bus(j,3)))+2*bus(startpos,2)^2*By(startpos,endpos);
        J2(i,startpos)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3))+By(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3)));
        L2(i,startpos)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3))-By(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3)))-2*bus(startpos,2)^2*line(i,7)+bus(startpos,2)^2*By(startpos,endpos);
        
    end
    if (endpos>nswing)
        H2(i,endpos-1)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3))-By(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3)));
        N2(i,endpos-1)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3))+By(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3)));
        J2(i,endpos-1)=bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3))-By(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3)));
        L2(i,endpos-1)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3))-By(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3)));
    elseif (startpos<nswing)
        H2(i,endpos)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3))-By(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3)));
        N2(i,endpos)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3))+By(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3)));
        J2(i,endpos)=bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3))-By(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3)));
        L2(i,endpos)=-bus(startpos,2)*bus(endpos,2)*(Gy(startpos,endpos)*sin(bus(startpos,3)-bus(endpos,3))-By(startpos,endpos)*cos(bus(startpos,3)-bus(endpos,3)));
    end
end

H2=-H2;
N2=-N2;
J2=-J2;
L2=-L2;
% zhilu guan xi de Juzhen G0
G0(nline,178)=0;
G0(1:nline,1:178)=H2;
G0((nline+1):(2*nline),1:178)=J2;
G0(1:nline,(178+1):(2*178))=N2;
G0((nline+1):(2*nline),(178+1):(2*178))=L2;
% zhi lu yu zhu ru gonglv de guanxi juzhen
T0=G0*S0;


% read random Pg info. 9th de lei ji liang
pdfgen=csvread('pdfgen1.csv');
ngen=length(pdfgen(:,1));
PgPx=zeros(ngen,9);
PgQx=zeros(ngen,9);
for i=1:ngen
    % pdfgen(i,4) column 4 is each generator capacity
    % pdfgen(i,5) column 5 is no. of generators
    % pdfgen(i,6) column 6 is the availability
    % K1=pdfgen(i,7);
    PgPx(i,1:9)=NcalGCum(pdfgen(i,4), pdfgen(i,5), pdfgen(i,6));
    % PgPx(K1,10)=pdfgen(i,7);
    PgQx(i,1:9)=NcalGCum(pdfgen(i,3)/pdfgen(i,5), pdfgen(i,5), pdfgen(i,6));
    % PgQx(K1,10)=pdfgen(i,7);
end

% read random Pl info. 9th
pdfload=csvread('pdfload1.csv');
nload=length(pdfload(:,1));
PlPx=zeros(nload,9);
PlQx=zeros(nload,9);
for i=1:nload
    muiP=abs(pdfload(i,4));
    sgmiP=pdfload(i,5)/100*muiP;
    muiQ=abs(pdfload(i,3));
    sgmiQ=pdfload(i,5)/100*muiQ;
    % load is negative
    % K2=pdfload(i,6);
    PlPx(i,1:9)=NcalLCum(muiP,sgmiP);
    % PlPx(K2,10)=pdfload(i,6);
    PlQx(i,1:9)=NcalLCum(muiQ,sgmiQ);
    % PlQx(K2,10)=pdfload(i,6); 
end
            
% zhu ru de fa dian ji he fuhe de da juzhen

Pgx=zeros(2*(ngen),9);
Pgx(1:ngen,:)=PgPx;
Pgx((ngen+1):(2*ngen),:)=PgQx;

Plx=zeros((2*nload),9);
Plx(1:nload,:)=PlPx;
Plx((nload+1):(2*nload),:)=PlQx;

% zhilu de gonglv de Cum 9th.
% zhilu you gong guan yu fa dian ji zu
PLineX_p=zeros(nline,9);
for i=1:nline
    for j=1:(2*ngen)
        if (j<=ngen)
           nbg=pdfgen(j,7);
           if nbg<nswing
               aij=T0(i,nbg);
           else
               aij=T0(i,nbg-1);
           end
        elseif j>ngen 
            if pdfgen((j-ngen),7)<nswing
                aij=T0(i,(nbg-1+178));
            else
                aij=T0(i,(nbg-1+178-1));
            end
        end
        if (aij~=0)
            PLineX_p(i,1)=PLineX_p(i,1)+aij*Pgx(j,1);
            PLineX_p(i,2)=PLineX_p(i,2)+(aij^2)*Pgx(j,2);
            PLineX_p(i,3)=PLineX_p(i,3)+(aij^3)*Pgx(j,3);
            PLineX_p(i,4)=PLineX_p(i,4)+(aij^4)*Pgx(j,4);
            PLineX_p(i,5)=PLineX_p(i,5)+(aij^5)*Pgx(j,5);
            PLineX_p(i,6)=PLineX_p(i,6)+(aij^6)*Pgx(j,6);
            PLineX_p(i,7)=PLineX_p(i,7)+(aij^7)*Pgx(j,7);
            PLineX_p(i,8)=PLineX_p(i,8)+(aij^8)*Pgx(j,8);
            PLineX_p(i,9)=PLineX_p(i,9)+(aij^9)*Pgx(j,9);
        end
    end
end
%  fu he 
for i=1:nline
    for j=1:((2*nload))
        if (j<=nload)
           nbl=pdfload(j,6);
           if nbl<nswing
               aij=-T0(i,nbl);
           else
               aij=-T0(i,nbl-1);
           end
        elseif j>nload 
            if pdfload((j-nload),6)<nswing
                aij=-T0(i,(nbl-1+178));
            else
                aij=-T0(i,(nbl-1+178-1));
            end
        end
        if (aij~=0)
            PLineX_p(i,1)=PLineX_p(i,1)+aij*Plx(j,1);
            PLineX_p(i,2)=PLineX_p(i,2)+(aij^2)*Plx(j,2);
            PLineX_p(i,3)=PLineX_p(i,3)+(aij^3)*Plx(j,3);
            PLineX_p(i,4)=PLineX_p(i,4)+(aij^4)*Plx(j,4);
            PLineX_p(i,5)=PLineX_p(i,5)+(aij^5)*Plx(j,5);
            PLineX_p(i,6)=PLineX_p(i,6)+(aij^6)*Plx(j,6);
            PLineX_p(i,7)=PLineX_p(i,7)+(aij^7)*Plx(j,7);
            PLineX_p(i,8)=PLineX_p(i,8)+(aij^8)*Plx(j,8);
            PLineX_p(i,9)=PLineX_p(i,9)+(aij^9)*Plx(j,9);
        end
    end
end

% zhilu wu gong guan yu fa dian ji zu
PLineX_q=zeros(nline,9);
for i=1:nline
    for j=1:(2*ngen)
        if (j<=ngen)
           nbg=pdfgen(j,7);
           if nbg<nswing
               aij=T0((i+178),nbg);
           else
               aij=T0((i+178),nbg-1);
           end
        elseif j>ngen 
            if pdfgen((j-ngen),7)<nswing
                aij=T0((i+178),(nbg-1+178));
            else
                aij=T0((i+178),(nbg-1+178-1));
            end
        end
        if (aij~=0)
            PLineX_q(i,1)=PLineX_q(i,1)+aij*Pgx(j,1);
            PLineX_q(i,2)=PLineX_q(i,2)+(aij^2)*Pgx(j,2);
            PLineX_q(i,3)=PLineX_q(i,3)+(aij^3)*Pgx(j,3);
            PLineX_q(i,4)=PLineX_q(i,4)+(aij^4)*Pgx(j,4);
            PLineX_q(i,5)=PLineX_q(i,5)+(aij^5)*Pgx(j,5);
            PLineX_q(i,6)=PLineX_q(i,6)+(aij^6)*Pgx(j,6);
            PLineX_q(i,7)=PLineX_q(i,7)+(aij^7)*Pgx(j,7);
            PLineX_q(i,8)=PLineX_q(i,8)+(aij^8)*Pgx(j,8);
            PLineX_q(i,9)=PLineX_q(i,9)+(aij^9)*Pgx(j,9);
        end
    end
end
%  fu he 
for i=1:nline
    for j=1:((2*nload))
        if (j<=nload)
           nbl=pdfload(j,6);
           if nbl<nswing
               aij=-T0((i+178),nbl);
           else
               aij=-T0((i+178),nbl-1);
           end
        elseif j>nload 
            if pdfload((j-nload),6)<nswing
                aij=-T0((i+178),(nbl-1+178));
            else
                aij=-T0((i+178),(nbl-1+178-1));
            end
        end
        if (aij~=0)
            PLineX_q(i,1)=PLineX_q(i,1)+aij*Plx(j,1);
            PLineX_q(i,2)=PLineX_q(i,2)+(aij^2)*Plx(j,2);
            PLineX_q(i,3)=PLineX_q(i,3)+(aij^3)*Plx(j,3);
            PLineX_q(i,4)=PLineX_q(i,4)+(aij^4)*Plx(j,4);
            PLineX_q(i,5)=PLineX_q(i,5)+(aij^5)*Plx(j,5);
            PLineX_q(i,6)=PLineX_q(i,6)+(aij^6)*Plx(j,6);
            PLineX_q(i,7)=PLineX_q(i,7)+(aij^7)*Plx(j,7);
            PLineX_q(i,8)=PLineX_q(i,8)+(aij^8)*Plx(j,8);
            PLineX_q(i,9)=PLineX_q(i,9)+(aij^9)*Plx(j,9);
        end
    end
end
% zhilu gonglv de Cum

PLineX((2*nline),9)=0;
PLineX(1:nline,:)=PLineX_p;
PLineX((nline+1):(2*nline),:)=PLineX_q;
% zhilu you gong de zhongxinju he Gram Charles index
PLineU_p=zeros(nline,9);
for i=1:nline
    PLineU_p(i,:)=CalU1(PLineX(i,1), PLineX(i,2), PLineX(i,3), PLineX(i,4), PLineX(i,5), PLineX(i,6), PLineX(i,7), PLineX(i,8), PLineX(i,9));
end

PLineC_p=zeros(nline,9);
for i=1:nline
    smg=sqrt(PLineX(i,2));
    PLineC_p(i,:)=CalC1(smg, PLineU_p(i,1), PLineU_p(i,2), PLineU_p(i,3), PLineU_p(i,4), PLineU_p(i,5), PLineU_p(i,6), PLineU_p(i,7), PLineU_p(i,8), PLineU_p(i,9));
end 

% zhilu wu gong de zhongxinju he Gram Charles index
PLineU_q=zeros(nline,9);
for i=1:nline
    PLineU_q(i,:)=CalU1(PLineX((i+nline),1), PLineX((i+nline),2), PLineX((i+nline),3), PLineX((i+nline),4), PLineX((i+nline),5), PLineX((i+nline),6), PLineX((i+nline),7), PLineX((i+nline),8), PLineX((i+nline),9));
end

PLineC_q=zeros(nline,9);
for i=1:nline
    smg=sqrt(PLineX(i,2));
    PLineC_q(i,:)=CalC1(smg, PLineU_q(i,1), PLineU_q(i,2), PLineU_q(i,3), PLineU_q(i,4), PLineU_q(i,5), PLineU_q(i,6), PLineU_q(i,7), PLineU_q(i,8), PLineU_q(i,9));
end 
% zhilu de  Gram Charles index
PLineC(2*nline,9)=0;
PLineC(1:nline,:)=PLineC_p;
PLineC((nline+1):(2*nline),:)=PLineC_q;

tt=cputime-t