pes2014-basic/test.m

%% Test
clear
arraymaxPDError=zeros(54,1);
for HH=0:53
    
    arraymaxPDError(HH+1)=-100;
    for CC=1:250
        tic
        [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]= ...
            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');
        %pf('ieee10471PG.dat');
        
        %PVi<EFBFBD><EFBFBD>ѹ<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        %PVu<EFBFBD><EFBFBD>ѹ<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
        Volt;
        UAngel*180/3.1415926;
        %Precision=Precision/1000000;
        %Precision=Precision*1000;
        %% ͨ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>PG
        AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
        PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
        
        %% <EFBFBD><EFBFBD>ֵ
        PG0=PG;
        PD0=PD;
        %%
        PG0(Balance)=PGBal(Balance);
        %%
        [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,HH);
        Gap=(Init_L*Init_Z'-Init_U*Init_W');
        KK=0;
        plotGap=zeros(1,50);
        ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
        kmax=600;
        %% 20120523 <EFBFBD><EFBFBD>ʱ
        QD_NON_ZERO=QD(PD==0 & QD~=0);
        QD_NON_ZERO_IND=find(PD==0 & QD~=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);
            %%
            L_1Z=diag(Init_Z./Init_L);
            U_1W=diag(Init_W./Init_U);
            %% <EFBFBD>γɺ<EFBFBD>ɭ<EFBFBD><EFBFBD>
            deltdeltF=func_deltdeltF(PVi,wG,wD,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);
            %%  <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>deltF
            deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,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,Loadi);
            Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND);
            Ly=Mat_H;
            Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK);
            Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK);
            Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
            %%YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
            %% <EFBFBD><EFBFBD>ʼ<EFBFBD>ⷽ<EFBFBD><EFBFBD>
            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]=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,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);
        ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi);
        %DrawGap(plotGap);
        %%
        %Volt=full(Volt');
        %PD=full(PD);
        %% ͳ<EFBFBD><EFBFBD>PD<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        maxPDError=max(abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  ));
        if maxPDError>arraymaxPDError(HH+1)
            arraymaxPDError(HH+1)=maxPDError;
        end
        toc
    end
end