随机测试缺负荷

Signed-off-by: facat <dmy@dmy-PC.(none)>

FormH.m

@@ -1,8 +1,8 @@
 function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND)
 %% 
-% QDcos=textread('1047glys.txt');
+QDcos=textread('D:\Project\×îС»¯³±Á÷\×îС³±Á÷ËãÀý\Ïɺ£919PDQDglys.txt');
-% QD(QD~=0)=PD(QD~=0)./tan(QDcos);
+QD(QD~=0)=PD(QD~=0)./tan(QDcos);
-% QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
+%QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
 %% 
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
 dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';

FormLw.asv

@@ -1,12 +1,14 @@
 function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK)
-
+%KK=999;
 PU=1*GenU(:,2);%发电机有功上界
 QU=PVQU(:,1);%发电机无功上界
-VoltU=(1.2+1/ekk)*ones(1,Busnum);
+%VoltU=(1.5+1/exp(KK))*ones(1,Busnum);
+VoltU=1.5*ones(1,Busnum);
 PDU=PD0(Loadi);
-PDU(PDU>0)=1.2*PDU(PDU>0);
+% PDU(PDU>0)=1200*PDU(PDU>0);
-PDU(PDU<0)=0.8*PDU(PDU<0);
+% PDU(PDU<0)=-800*PDU(PDU<0);
-PDU(PDU==0)=.2;
+% PDU(PDU==0)=200;
+PDU=1000000*ones(length(Loadi),1);
 t1=([PU',QU',PDU',VoltU])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;

FormLw.m

@@ -1,12 +1,14 @@
 function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK)
-%KK=999;
+KK=999;
-PU=1*GenU(:,2);%发电机有功上界
+PU=GenU(:,2);%发电机有功上界
 QU=PVQU(:,1);%发电机无功上界
 VoltU=(1.2+1/exp(KK))*ones(1,Busnum);
+%VoltU=10*ones(1,Busnum);
 PDU=PD0(Loadi);
-PDU(PDU>0)=1.20*PDU(PDU>0);
+PDU(PDU>0)=1.200*PDU(PDU>0);
-PDU(PDU<0)=.80*PDU(PDU<0);
+PDU(PDU<0)=0.800*PDU(PDU<0);
-PDU(PDU==0)=.20;
+PDU(PDU==0)=0.400;
+%PDU=10*ones(length(Loadi),1);
 t1=([PU',QU',PDU',VoltU])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;

FormLz.m

@@ -1,12 +1,14 @@
 function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK)
-%KK=999;
+KK=999;
-PL=1*GenL(:,2);%发电机有功下界
+PL=GenL(:,2);%发电机有功下界
 QL=PVQL(:,1);%发电机无功下界
 VoltL=(0.8-1/exp(KK))*ones(1,Busnum);
+%VoltL=-10*ones(1,Busnum);
 PDL=PD0(Loadi);
-PDL(PDL>0)=.80*PDL(PDL>0);
+PDL(PDL>0)=0.800*PDL(PDL>0);
-PDL(PDL<0)=1.20*PDL(PDL<0);
+PDL(PDL<0)=1.200*PDL(PDL<0);
-PDL(PDL==0)=-.20;
+PDL(PDL==0)=-0.400;
+%PDL=-10*ones(length(Loadi),1);
 t1=([PL',QL',PDL',VoltL])';
 t2=Mat_G-Init_L'-t1;
 Lz=t2;

FormYY.m

@@ -1,4 +1,4 @@
-function YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx)
+function YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx)
 YY=[
   Lx;
   -Ly;

OPF.asv

@@ -1,10 +1,18 @@
 tic
 clear
-[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(c:/iee.txt);
+[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\最小化潮流\最小潮流算例\金湖924（2-1）_0.5_85%.txt'); 
+%pf('c:/file31.txt'); 
+%pf('ieee10471PG.dat');
+
 %PVi电压节点序号
 %PVu电压节点电压标幺值
 Volt;
 UAngel*180/3.1415926;
+%% 通过潮流计算PG
+AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
+dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
+PG0()
 %% 初值
 PG0=PG;
 PD0=PD;
@@ -13,7 +21,7 @@ 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;
+kmax=60;
 %% 20120523 临时
 QD_NON_ZERO=QD(PD==0 & QD~=0);
 QD_NON_ZERO_IND=find(PD==0 & QD~=0);
@@ -48,8 +56,8 @@ while(abs(Gap)>Precision)
     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);
+    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi);
+    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);
     %% 开始解方程
@@ -64,7 +72,9 @@ fprintf('
 ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
 DrawGap(plotGap);
 %% 
-% Volt=full(Volt');
+Volt=full(Volt');
-% PD=full(PD);
+PD=full(PD);
+%% 统计PD误差
+abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  )
 toc  
 

OPF.m

@@ -1,13 +1,24 @@
 tic
 clear
-[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('e:/算例/standard.txt');
+[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]= ...
-%PVi电压节点序号
+pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); 
-%PVu电压节点电压标幺值
+%pf('c:/file31.txt'); 
+%pf('ieee10471PG.dat');
+
+%% 计算功率因数
+atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0))
 Volt;
 UAngel*180/3.1415926;
+%% 通过潮流计算PG
+AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
+PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
+
 %% ³õÖµ
 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);
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
@@ -64,7 +75,9 @@ fprintf('
 ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
 DrawGap(plotGap);
 %% 
-% Volt=full(Volt');
+%Volt=full(Volt');
-% PD=full(PD);
+%PD=full(PD);
+%% 统计PD误差
+maxPDError=max(abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  ));
 toc  
 

OPF_Init.m

@@ -1,4 +1,4 @@
-function [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)
+function [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)
 Loadi=find(QD~=0);
 %Loadi=[1:Busnum]';
 RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1; %约束条件数
@@ -17,13 +17,13 @@ tPL=sparse(GenL(:,2));%
 tQL=sparse(PVQL(:,1));% 无功下限
 PG(PGi)=(tPU+tPL)/2;
 QG(PVi)=(tQU+tQL)/2;
-wG=ones(size(PGi,1),1);
+wG=0.1*ones(size(PGi,1),1);
 randInt=randperm(size(Loadi,1));
-randPDind=randInt(1:0);
+randPDind=randInt(1:HH);
-
+wD=0.1*ones(size(Loadi,1),1);
-wD=ones(size(Loadi,1),1);
 wD(randPDind)=0;%一些负荷不约束
-%wD(12)=0;
+%wD(7)=0;
+% wD(11)=0;
 PD=1*PD0;
 
 end

SolveIt.m

@@ -1,7 +1,8 @@
 function 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)
 LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx);
-H=-deltdeltF+ddh;%+ddg*(Init_Z'+Init_W');
+H=-deltdeltF+ddh;
-t1=diag(Init_L.\Init_Z-Init_U.\Init_W);
+%t1=diag(Init_L.\Init_Z-Init_U.\Init_W);
+t1=diag(Init_Z./Init_L-Init_W./Init_U);
 t2=-deltG*( t1 )*deltG';
 aa=[
     (H+t2),deltH;

func_ddh.m

@@ -1,4 +1,4 @@
-function ddh=func_ddh3(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount)
+function ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount)
 %决定用循环重写
 %ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
 
@@ -42,8 +42,11 @@ ddQdTdV=t1+t2+t3+t4;
 t1=Y.*sin(mat_INV_AngleIJ)*diag(yQ);
 t2=diag(yQ)*Y.*sin(mat_AngleIJ);
 ddQdVdV=t1+t2;
-t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV ;
+% t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV;
-   ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT;
+%    ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT;
+% ];
+t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
+   ddPdTdV+ddQdTdV,ddPdTdT+ddQdTdT;
 ];
 sizePGi=size(PGi,1);
 sizePVi=size(PVi,1);

func_deltF.m

@@ -4,7 +4,7 @@ t1=2*wG.*(PG(PGi) - PG0(PGi) );
 t2=2*wD.*(PD(Loadi)-PD0(Loadi));
 deltF=[
     sparse(t1);
-    sparse(size(PVi,1),size(PVi,2));
+    sparse(size(PVi,1),1);
     sparse(t2);
     sparse(2*Busnum,1);
 ];

test.m

@@ -0,0 +1,97 @@
+%% 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\最小化潮流\最小潮流算例\仙海919.txt');
+        %pf('c:/file31.txt');
+        %pf('ieee10471PG.dat');
+        
+        %PVi电压节点序号
+        %PVu电压节点电压标幺值
+        Volt;
+        UAngel*180/3.1415926;
+        %Precision=Precision/1000000;
+        %Precision=Precision*1000;
+        %% 通过潮流计算PG
+        AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
+        PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
+        
+        %% 初值
+        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 临时
+        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;
+            %% 开始计算OPF
+            %% 形成等式约束的雅克比
+            deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
+            %% 形成不等式约束的雅克比
+            deltG=func_deltG(Busnum,PVi,PGi,Loadi);
+            %%
+            L_1Z=diag(Init_Z./Init_L);
+            U_1W=diag(Init_W./Init_U);
+            %% 形成海森阵
+            deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount);
+            %% 形成ddHy
+            ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
+            %%  开始构建ddg
+            ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi);
+            %%  开始构建deltF
+            deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi);
+            
+            %% 形成方程矩阵
+            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);
+            %% 开始解方程
+            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);
+            %%取各分量
+            [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('迭代次数%d\n',KK);
+        ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi);
+        %DrawGap(plotGap);
+        %%
+        %Volt=full(Volt');
+        %PD=full(PD);
+        %% 统计PD误差
+        maxPDError=max(abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  ));
+        if maxPDError>arraymaxPDError(HH+1)
+            arraymaxPDError(HH+1)=maxPDError;
+        end
+        toc
+    end
+end
+
+