nce3.m

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

clear
clear global
% global variables
glo_var
% Edgeworth's Expansion Theory

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. 9th
pdfgen=csvread('pdfgen1.csv');
ngen=length(pdfgen(:,1));
Pgx=zeros(ngen,3);
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,:)=cal3GCum(pdfgen(i,4), pdfgen(i,5), pdfgen(i,6));
end

% read random Pl info. 9th
pdfload=csvread('pdfload.csv');
nload=length(pdfload(:,1));
Plx=zeros(nload,3);
for i=1:nload
    mui=pdfload(i,4);
    sgmi=pdfload(i,5)/100*mui;
    % load is negative
    Plx(i,:)=cal3LCum(mui,sgmi);
end

Pglx=zeros(ngen+nload,3);
Pglx(1:ngen,:)=Pgx;
Pglx(ngen+1:ngen+nload,:)=Plx;

PLineX=zeros(nline,3);
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);
        end
    end
end

PLineU=zeros(nline,3);
for i=1:nline
    PLineU(i,:)=CalU3(PLineX(i,1), PLineX(i,2), PLineX(i,3));
end

PLineC=zeros(nline,3);
for i=1:nline
    smg=sqrt(PLineX(i,2));
    PLineC(i,:)=CalC3(smg, PLineU(i,1), PLineU(i,2), PLineU(i,3));
end 
        
tt=cputime-t