mcwscc2.m

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

clear
clear global
% global variables
glo_var

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);

% Mpd is the relation matrix between Pline and Delta(voltage angle)
Mpd=formMpd(bus, line, nswing);
% relationship matrix between Pline and Pbus
Mind=Mpd*inv(Yred);

% read random Pg info.
pdfgen=csvread('pdfgen1.csv');
ngen=length(pdfgen(:,1));
Pgx=zeros(ngen,6);
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
    Pgx(i,:)=calGCum(pdfgen(i,4), pdfgen(i,5), pdfgen(i,6));
end

% read random new generation Pgn info.
pdfngen=csvread('pdfnew.csv');
nnewgen=length(pdfngen(:,1));
Pgnx=zeros(nnewgen,6);
for i=1:nnewgen
    % pdfngen(i,2) column 2 is each generator capacity
    % pdfngen(i,3) column 3 is no. of generators
    % pdfngen(i,4) column 4 is the availability
    Pgnx(i,:)=calGCum(pdfngen(i,2), pdfngen(i,3), pdfngen(i,4));
end
    
% read random Pl info.
pdfload=csvread('pdf41load.csv');
nload=length(pdfload(:,1));
for i=1:nload
    mui=pdfload(i,4);
    sgmi=pdfload(i,5)/100*mui;
    % load is negative
    Plx(i,:)=calLCum(mui,sgmi);
end

Pglx=zeros(ngen+nload,6);
Pglx(1:ngen,:)=Pgx;
Pglx(ngen+1:ngen+nnewgen,:)=Pgnx;
Pglx(ngen+nnewgen+1:ngen+nnewgen+nload,:)=Plx;

PLineX=zeros(nline,6);
for i=1:nline
    for j=1:ngen+nnewgen+nload
        if (j<=ngen)
            nbg=pdfgen(j,7);
            aij=findInd(i, nbg, Mind, nswing);
        elseif (j>ngen)&(j<=ngen+nnewgen) % (new generation)
            nbng=pdfngen(j-ngen,5);
            aij=findInd(i, nbng, Mind, nswing);
        else 
            nbl=pdfload(j-ngen-nnewgen,6);
            % load is negative
            aij=-findInd(i, nbl, Mind, nswing);
        end
        if (aij~=0)
            PLineX(i,1)=PLineX(i,1)+aij*Pglx(j,1);
            PLineX(i,2)=PLineX(i,2)+aij^2*Pglx(j,2);
            PLineX(i,3)=PLineX(i,3)+aij^3*Pglx(j,3);
            PLineX(i,4)=PLineX(i,4)+aij^4*Pglx(j,4);
            PLineX(i,5)=PLineX(i,5)+aij^5*Pglx(j,5);
            PLineX(i,6)=PLineX(i,6)+aij^6*Pglx(j,6);
        end
    end
end

PLineC=zeros(nline,6);
for i=1:nline
    PLineC(i,:)=CalC(PLineX(i,1), PLineX(i,2), PLineX(i,3), PLineX(i,4), PLineX(i,5), PLineX(i,6));
end
        
tt=cputime-t