pes2014/OPF.m

tic
clc
clear
thesis=ForThesis(4,8);
for II=1:4
    [kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL,Liner,Linex,Lineb,Transforr,Transforx,Transfork0]= ...
        pf('c:/newFIle.txt');
    %pf('D:\Project\<EFBFBD><EFBFBD><EFBFBD>㽵<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ŀ\<EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\<EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\ԭʼ\<EFBFBD><EFBFBD>ͷվ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>922<EFBFBD><EFBFBD>3-1<EFBFBD><EFBFBD>_0.5_120%.txt');
    %pf('D:\Project\<EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\<EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\<EFBFBD>ɺ<EFBFBD>919.txt');
    %pf('c:/file31.txt');
    
    
    %% <EFBFBD><EFBFBD><EFBFBD>㹦<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    Loadi=QD~=0 | PD~=0;
    PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2))
    %%
    Volt;
    UAngel*180/3.1415926;
    %% ͨ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>PG
    AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
    PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
    QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';
    %% <EFBFBD><EFBFBD>ֵ-<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
    PG0=PG;
    QG0=QG;
    PD0=PD;
    QD0=QD;
    PDReal=PD;%<EFBFBD><EFBFBD>ֵ
    QDReal=QD;%<EFBFBD><EFBFBD>ֵ
    %PD0(12)=PD0(12)+0.001;
    %%
    PG0(Balance)=PGBal(Balance);
    PG(Balance)=PGBal(Balance);
    QG0(Balance)=QGBal(Balance);
    QG0(PVi)=QGBal(PVi);
    QG(PVi)=QGBal(PVi);
    %%
    PVi=zeros(0,0);
    [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD);
    Gap=(Init_L*Init_Z'-Init_U*Init_W');
    KK=0;
    plotGap=zeros(1,60);
    ContrlCount=size(PGi,1)+size(Loadi,1)*2+Busnum*2;
    kmax=60;
    %% 20120523 <EFBFBD><EFBFBD>ʱ
    QD_NON_ZERO=QD(PD==0 & QD~=0);
    QD_NON_ZERO_IND=find(PD==0 & QD~=0);
    %%
    Precision=Precision/10;
    
    %% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
    QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
%     PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
%     QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
    %% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    %[noDataTransNum noDataTransCapacity noDataTransPowerFactor]=ReadNoDataTrans('C:/b/<EFBFBD><EFBFBD><EFBFBD><EFBFBD>911_2751267_2012-09-05/iPso_<EFBFBD><EFBFBD><EFBFBD><EFBFBD>911_2751267_2012-09-05_<EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD>޸<EFBFBD><EFBFBD><EFBFBD>.txt');
    noDataTransCapacity=0;
    
    while(abs(Gap)>Precision)
        if KK>kmax
            break;
        end
        plotGap(KK+1)=Gap;
        Init_u=Gap/2/RestraintCount*CenterA;
        AngleIJMat=0;
        %% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>OPF
        %% <EFBFBD>γɵ<EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ſ˱<EFBFBD>
        deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
        %% <EFBFBD>γɲ<EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ſ˱<EFBFBD>
        deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD);
        %%
        L_1Z=diag(Init_Z./Init_L);
        U_1W=diag(Init_W./Init_U);
        %% <EFBFBD>γɺ<EFBFBD>ɭ<EFBFBD><EFBFBD>
        deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
        %% <EFBFBD>γ<EFBFBD>ddHy
        ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
        %%  <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ddg
        ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD);
        %%  <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>deltF
        deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
        
        %% <EFBFBD>γɷ<EFBFBD><EFBFBD>̾<EFBFBD><EFBFBD><EFBFBD>
        Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
        Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
        Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi);
        Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND,Loadi);
        Ly=Mat_H;
        Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF);
        Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity);
        Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
        YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
        %% <EFBFBD><EFBFBD>ʼ<EFBFBD>ⷽ<EFBFBD><EFBFBD>
        fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> %d Gap %f\n',KK+1,plotGap(KK+1));
        XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
        %%ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
        [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi);
        Gap=(Init_L*Init_Z'-Init_U*Init_W');
        KK=KK+1;
    end
    fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>%d\n',KK);
    fprintf('Ŀ<EFBFBD><EFBFBD>ֵ%f\n',full(ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi)));
    % DrawGap(plotGap);
    %%
    %Volt=full(Volt');
    %PD=full(PD);
    %% ͳ<EFBFBD><EFBFBD>PD<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    % absPDLoad=abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  );
    absPDLoad=abs(  (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi)  );
    maxPDError=max(absPDLoad(absPDLoad<10))
    absQDLoad=abs(  (QD(Loadi)-QDReal(Loadi))./QDReal(Loadi)  );
    maxQDError=max(absQDLoad(absQDLoad<10))
    disp('index');
    %Loadi(absPDLoad==maxPDError);
    %% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    totalLoss=(sum(PG)-sum(PD(Loadi)))*100;
    fprintf('<EFBFBD>ܵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ%f(MW <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ)\n',full(totalLoss));
    fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ %f\n',full(totalLoss/sum(PG)));
    %% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    %Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel);
    thesis=thesis.AddPDQDPGQG(PD(Loadi),QD(Loadi),PG(Balance),QG(PVi));
    
end
PD(Loadi)=thesis.MeanPD();
QD(Loadi)=thesis.MeanQD();
PG(Balance)=thesis.MeanPG();
QG(PVi)=thesis.MeanQG();
thesis.MaxDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi));
thesis.StatDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))
toc