mcwscc1.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 Pl info.
pdfload=csvread('pdfload.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+nload,:)=Plx;

PLineX=zeros(nline,6);
for i=1:nline
    for j=1:ngen+nload
        if (j<=ngen)
            nbg=pdfgen(j,7);
            aij=findInd(i, nbg, Mind, nswing);
        else % (j>=ngen)
            nbl=pdfload(j-ngen,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