From e515e5fd4a25333433587f22840c36c84c430cf7 Mon Sep 17 00:00:00 2001 From: facat Date: Fri, 25 May 2012 17:16:59 +0800 Subject: [PATCH] =?UTF-8?q?=E5=88=A0=E6=8E=89=E5=A4=9A=E4=BD=99=E7=9A=84?= =?UTF-8?q?=E5=87=BD=E6=95=B0=E5=8F=82=E6=95=B0=EF=BC=8C=E7=BB=9F=E4=B8=80?= =?UTF-8?q?=E4=BA=86ContrlCount?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- FormG.m | 14 +------------- FormH.m | 2 +- FormLw.m | 2 +- FormLz.m | 2 +- OPF.asv | 34 ++++++++++++++++++---------------- OPF.m | 14 +++++++------- OPF_Init.m | 8 +------- func_ddh3.m | 4 ++-- func_deltF.m | 11 ++--------- func_deltH.asv | 24 ++++++------------------ func_deltH.m | 2 +- func_deltdeltF.asv | 14 ++++++++------ func_deltdeltF.m | 4 ++-- 13 files changed, 51 insertions(+), 84 deletions(-) diff --git a/FormG.m b/FormG.m index e336b9a..4f8b4ea 100644 --- a/FormG.m +++ b/FormG.m @@ -1,21 +1,9 @@ function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,Loadi) -%t1=PG(PVi); -%GP=t1;%发电机P -%GP=[4.5 4.5]'; -%%线路 -%发电机Q -% t1=Volt'*Volt; -% t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -% t3=t1.*t2; -% t4=sum(t3,2);%发电机Q -%GQ=t4; + Mat_G=[ - %GP; PG(PGi); QG(PVi); sparse(PD(Loadi)); - %GQ(PVi); - %[0 1.45]'; Volt'; ]; end \ No newline at end of file diff --git a/FormH.m b/FormH.m index e0f68a6..87a5378 100644 --- a/FormH.m +++ b/FormH.m @@ -1,4 +1,4 @@ -function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND) +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'); QD(QD~=0)=PD(QD~=0)./tan(QDcos); diff --git a/FormLw.m b/FormLw.m index c2312dd..2029450 100644 --- a/FormLw.m +++ b/FormLw.m @@ -1,4 +1,4 @@ -function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi) +function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi) PU=1*GenU(:,2);%发电机有功上界 QU=PVQU(:,1);%发电机无功上界 diff --git a/FormLz.m b/FormLz.m index b85ce45..31263fb 100644 --- a/FormLz.m +++ b/FormLz.m @@ -1,4 +1,4 @@ -function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi) +function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi) PL=1*GenL(:,2);%发电机有功下界 QL=PVQL(:,1);%发电机无功下界 diff --git a/OPF.asv b/OPF.asv index 73b28d1..9e397ab 100644 --- a/OPF.asv +++ b/OPF.asv @@ -1,6 +1,6 @@ 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('ieee3001PG.dat'); +[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('ieee10471PG.dat'); %PVi电压节点序号 %PVu电压节点电压标幺值 Volt; @@ -8,11 +8,11 @@ UAngel*180/3.1415926; %% 初值 PG0=PG; PD0=PD; -[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,randPDind]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0); +[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; plotGap=zeros(1,50); -ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; +ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2; kmax=600; %% 20120523 临时 QD_NON_ZERO=QD(PD==0 & QD~=0); @@ -27,42 +27,44 @@ while(abs(Gap)>Precision) AngleIJMat=0; %% 开始计算OPF %% 形成等式约束的雅克比 - deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle); + deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi); %% 形成不等式约束的雅克比 - deltG=func_deltG(Busnum,PVi,PGi); + deltG=func_deltG(Busnum,PVi,PGi,Loadi); %% L_1Z=diag(Init_Z./Init_L); U_1W=diag(Init_W./Init_U); %% 形成海森阵 - deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0); + deltdeltF=func_deltdeltF(Busnum,PVi,PGi,wG,wD,ContrlCount); %% 形成ddHy - ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle); + ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi); %% 开始构建ddg - ddg=func_ddg(PGi,PVi,Busnum,RestraintCount); + ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi); %% 开始构建deltF - deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum); + deltF=func_deltF(PG,PVi,GenC,ContrlCount,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); + Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,Loadi); Mat_H=FormH(Busnum,GB,AngleIJMat,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); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0); + Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi); + Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi); Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); YY=FormYY1(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); + 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]=AssignXX1(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); + [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) +ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi) DrawGap(plotGap); -Volt=Volt'; +%% +Volt=full(Volt'); +PD=full(PD); toc diff --git a/OPF.m b/OPF.m index dc3c6f8..e031e8e 100644 --- a/OPF.m +++ b/OPF.m @@ -27,29 +27,29 @@ while(abs(Gap)>Precision) AngleIJMat=0; %% 开始计算OPF %% 形成等式约束的雅克比 - deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle,Loadi); + 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(Busnum,GenC,PVi,PGi,wG,wD,PD0); + deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount); %% 形成ddHy - ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi); + ddh=func_ddh3(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,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi); + 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,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND); + 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,randPDind,Loadi); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi); + Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi); + Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi); Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx); %% 开始解方程 diff --git a/OPF_Init.m b/OPF_Init.m index 7f18f70..f0c4745 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -19,16 +19,10 @@ PG(PGi)=(tPU+tPL)/2; QG(PVi)=(tQU+tQL)/2; wG=ones(size(PGi,1),1); randInt=randperm(size(Loadi,1)); -%asc_randInt=sort(randInt,2,'ascend' ); -%randPDind=randInt(1:200); - -randPDind=randInt(1:0); +randPDind=randInt(1:262); wD=ones(size(Loadi,1),1); wD(randPDind)=0;%一些负荷不约束 -%wD(Balance)=0; -%PD0(randPDind)=0; PD=1*PD0; -%PD(PD==0)=.2; end \ No newline at end of file diff --git a/func_ddh3.m b/func_ddh3.m index 0afb739..22ed8fb 100644 --- a/func_ddh3.m +++ b/func_ddh3.m @@ -1,6 +1,6 @@ -function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi) +function ddh=func_ddh3(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; +%ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2; %AngleIJ=AngleIJMat-angle(GB); mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); diff --git a/func_deltF.m b/func_deltF.m index db6a5a2..3cf4c06 100644 --- a/func_deltF.m +++ b/func_deltF.m @@ -1,12 +1,5 @@ -function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi) -%t1=PG(setdiff(PVi,Balance)); -% t2=Volt'*Volt; -% t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -% t4=t2.*t3; -% t5=sum(t4,2); -% PBal=t5(Balance); -% PPG=([PQ(1),PBal])';%暂时用土办法处理一下 -%% +function deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi) + t1=2*wG.*(PG(PGi) - PG0(PGi) ); t2=2*wD.*(PD(Loadi)-PD0(Loadi)); deltF=[ diff --git a/func_deltH.asv b/func_deltH.asv index 56dfec9..317b2f3 100644 --- a/func_deltH.asv +++ b/func_deltH.asv @@ -1,20 +1,8 @@ -function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi) -dH_dPg=zeros(size(PGi,1)+size(PVi,1),2*Busnum); +function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi) -for I=1:size(PGi,1) - %dH_dPg(I,2*PVi(I)-1)=-1;王锡凡书上的公式 - %dH_dPg(I,2*PGi(I)-1)=1;暂时改一下20111227 - dH_dPg(I,PGi(I))=1; -end - -dH_dQr=zeros(size(PVi,1)+size(PGi,1),2*Busnum); -for I=1:size(PVi,1) - %dH_dQr(I,2*PVi(I))=-1;王锡凡书上的公式 - %dH_dQr(I,2*PVi(I))=1;暂时改一下20111227 - dH_dQr(I,PVi(I))=1; -end -Angle=angle(GB); -dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat); %形成雅克比矩阵 -%deltH=[dH_dPg;dH_dQr;dH_dx'];%dH_dx 需要使用一下转置 暂时改一下 -deltH=[dH_dPg;dH_dQr;dH_dx']; +dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum); +dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,1),2*Busnum); +dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)]; +dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c); %形成雅克比矩阵 +deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx']; end \ No newline at end of file diff --git a/func_deltH.m b/func_deltH.m index 871346d..317b2f3 100644 --- a/func_deltH.m +++ b/func_deltH.m @@ -1,4 +1,4 @@ -function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle,Loadi) +function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi) dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum); dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,1),2*Busnum); diff --git a/func_deltdeltF.asv b/func_deltdeltF.asv index 3a69884..87f3a1a 100644 --- a/func_deltdeltF.asv +++ b/func_deltdeltF.asv @@ -1,9 +1,11 @@ -function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0) +function deltdeltF=func_deltdeltF(Busnum,PVi,PGi,wG,wD,ContrlCount) -ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; %P,Q,Volt theta这些控制变量数 +ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数 C=[wG' zeros(size(PVi))' wD']; -% deltdeltF=[diag(C)*2,zeros(size(C,2),ContrlCount-size(C,2)); -% zeros(ContrlCount-size(C,2),ContrlCount); -% ]; - +sizeC=size(C,2); +diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC); +deltdeltF=[ + diagC*2,sparse(sizeC,ContrlCount-sizeC); + sparse(ContrlCount-sizeC,ContrlCount); + ]; end \ No newline at end of file diff --git a/func_deltdeltF.m b/func_deltdeltF.m index db1db61..a4317a3 100644 --- a/func_deltdeltF.m +++ b/func_deltdeltF.m @@ -1,6 +1,6 @@ -function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0) +function deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount) -ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数 +%ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数 C=[wG' zeros(size(PVi))' wD']; sizeC=size(C,2); diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC);