From b071a58e33d47c44908bf9c9f300641500c6c2e6 Mon Sep 17 00:00:00 2001 From: dmy Date: Wed, 14 May 2014 16:37:56 +0800 Subject: [PATCH] =?UTF-8?q?=E5=AE=8C=E6=88=90=E4=BA=86QD=EF=BC=8C=E5=B9=B6?= =?UTF-8?q?=E4=BF=AE=E5=A4=8D=E4=BA=86=E5=AE=8C=E6=88=90PD=E6=97=B6?= =?UTF-8?q?=E7=9A=84=E4=B8=80=E4=B8=AA=E6=B1=82=E5=AF=BC=E9=94=99=E8=AF=AF?= =?UTF-8?q?=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dmy --- FormG.m | 10 ++++++---- FormLw.m | 5 +++-- FormLz.m | 6 ++++-- Modification.m | 5 +++-- OPF.m | 18 +++++++++++------- OPF_Init.m | 2 +- func_ddh.m | 4 ++-- func_deltF.m | 1 + func_deltG.m | 33 ++++++++++++++++++++++++++++----- func_deltH.m | 3 ++- 10 files changed, 61 insertions(+), 26 deletions(-) diff --git a/FormG.m b/FormG.m index d8648ba..4f28c50 100644 --- a/FormG.m +++ b/FormG.m @@ -1,12 +1,14 @@ -function Mat_G=FormG(Volt,PD,QD,Loadi,Vbi,mVolt,RealVolt,bigM,sigma,PDbi,RealPD,RealQD) +function Mat_G=FormG(Volt,PD,QD,Loadi,Vbi,mVolt,RealVolt,bigM,sigma,PDbi,QDbi,RealPD,RealQD) Mat_G=[ - sparse(PD(Loadi))-bigM*PDbi-PD(Loadi)*sigma; - sparse(PD(Loadi))+bigM*PDbi+PD(Loadi)*sigma; - sparse(QD(Loadi)); + sparse(PD(Loadi))-bigM*PDbi-RealPD(Loadi)*sigma; + sparse(PD(Loadi))+bigM*PDbi+RealPD(Loadi)*sigma; + sparse(QD(Loadi))-bigM*QDbi-RealQD(Loadi)*sigma; + sparse(QD(Loadi))+bigM*QDbi+RealQD(Loadi)*sigma; Volt'-bigM*Vbi-mVolt'-sigma*RealVolt'; Volt'+bigM*Vbi+mVolt'+sigma*RealVolt'; Vbi; PDbi; + QDbi; ]; end \ No newline at end of file diff --git a/FormLw.m b/FormLw.m index 06b8737..c93fe39 100644 --- a/FormLw.m +++ b/FormLw.m @@ -1,4 +1,4 @@ -function Lw=FormLw(Mat_G,Init_U,Busnum,PD0,QD0,Loadi) +function [Lw,UpperLimit]=FormLw(Mat_G,Init_U,Busnum,PD0,QD0,Loadi,bigM) KK=999; VoltU=(1.1)*ones(1,Busnum); PDU=PD0(Loadi); @@ -17,7 +17,8 @@ realQD=QD0(Loadi); indQD=find(realQD>0); QDU(indQD(3:12:end))=1.55*realQD(indQD(3:12:end)); QDU(indQD(9:12:end))=1.05*realQD(indQD(9:12:end)); -t1=([0*PDU',1*ones(1,length(Loadi)),QDU',0*VoltU,5*ones(1,Busnum),1*ones(1,Busnum),1*ones(1,length(Loadi))])'; +t1=([0*PDU',bigM*ones(1,length(Loadi)),0*QDU',bigM*ones(1,length(Loadi)),0*VoltU,bigM*ones(1,Busnum),1*ones(1,Busnum),1*ones(1,length(Loadi)*2)])'; +UpperLimit=t1; t2=Mat_G+Init_U'-t1; Lw=t2; diff --git a/FormLz.m b/FormLz.m index 24f22ec..ff6a144 100644 --- a/FormLz.m +++ b/FormLz.m @@ -1,4 +1,4 @@ -function Lz=FormLz(Mat_G,Init_L,Busnum,PD0,QD0,Loadi) +function [Lz,LoweLimit]=FormLz(Mat_G,Init_L,Busnum,PD0,QD0,Loadi,bigM) VoltL=(0.9)*ones(1,Busnum); PDL=PD0(Loadi); PDL(PDL>0)=0.800*PDL(PDL>0); @@ -16,7 +16,9 @@ realQD=QD0(Loadi); indQD=find(realQD>0); QDL(indQD(3:12:end))=0.95*realQD(indQD(3:12:end)); QDL(indQD(9:12:end))=0.95*realQD(indQD(9:12:end)); -t1=([-1*ones(1,length(Loadi)),0*PDL',QDL',-5*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum),0*ones(1,length(Loadi))])'; +%t1=([-1*ones(1,length(Loadi)),0*PDL',-1*ones(1,length(Loadi)),0*QDL',-5*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum),0*ones(1,length(Loadi)*2)])'; +t1=([-bigM*ones(1,length(Loadi)),0*PDL',-bigM*ones(1,length(Loadi)),0*QDL',-bigM*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum),0*ones(1,length(Loadi)*2)])'; +LoweLimit=t1; t2=Mat_G-Init_L'-t1; Lz=t2; end \ No newline at end of file diff --git a/Modification.m b/Modification.m index 9b9fecc..c0b8e92 100644 --- a/Modification.m +++ b/Modification.m @@ -1,4 +1,4 @@ -function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi]=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,Vbi,PDbi) +function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=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,Vbi,PDbi,QDbi) AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU); %fprintf('AlphaP %f\n',full(AlphaP)); AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW); @@ -25,5 +25,6 @@ Volt=Volt+AlphaP*t(1:Busnum); Volt(Balance)=balVolt; UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum); Vbi=Vbi+AlphaP*t(2*Busnum+1:2*Busnum+Busnum)'; -PDbi=PDbi+AlphaP*t(2*Busnum+Busnum+1:end)'; +PDbi=PDbi+AlphaP*t(2*Busnum+Busnum+1:2*Busnum+Busnum+length(Loadi))'; +QDbi=QDbi+AlphaP*t(2*Busnum+Busnum+length(Loadi)+1:end)'; end \ No newline at end of file diff --git a/OPF.m b/OPF.m index a5aba6b..c075395 100644 --- a/OPF.m +++ b/OPF.m @@ -44,7 +44,7 @@ RealQD=QD0; Gap=(Init_L*Init_Z'-Init_U*Init_W'); KK=0; plotGap=zeros(1,60); -ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi); +ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2; kmax=600; Precision=Precision/1; sigma=0.03; @@ -57,7 +57,8 @@ RealVolt=Volt0; mVolt=Volt0.*(1+normrnd(0,sigma,1,length(Volt0))); Vbi=sparse(1*ones(Busnum,1)); PDbi=sparse(1*ones(length(Loadi),1)); -bigM=1; +QDbi=sparse(1*ones(length(Loadi),1)); +bigM=5; while(abs(Gap)>Precision) if KK>kmax break; @@ -84,23 +85,26 @@ while(abs(Gap)>Precision) %% 形成方程矩阵 Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1); Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1); - Mat_G=FormG(Volt,PD,QD,Loadi,Vbi,mVolt,RealVolt,bigM,sigma,PDbi,RealPD,RealQD); + Mat_G=FormG(Volt,PD,QD,Loadi,Vbi,mVolt,RealVolt,bigM,sigma,PDbi,QDbi,RealPD,RealQD); Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); Ly=Mat_H; - Lz=FormLz(Mat_G,Init_L,Busnum,PD0,QD0,Loadi); - Lw=FormLw(Mat_G,Init_U,Busnum,PD0,QD0,Loadi); + [Lz,LowerLimit]=FormLz(Mat_G,Init_L,Busnum,PD0,QD0,Loadi,bigM);%下界 + [Lw,UpperLimit]=FormLw(Mat_G,Init_U,Busnum,PD0,QD0,Loadi,bigM);%上界 Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); + if any(LowerLimit>UpperLimit) + warning('不等式设置错误'); + end YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); %% 开始解方程 fprintf('迭代次数 %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); %%取各分量 [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,Vbi,PDbi]=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,Vbi,PDbi); + [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=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,Vbi,PDbi,QDbi); Gap=(Init_L*Init_Z'-Init_U*Init_W'); KK=KK+1; end fprintf('目标函数: %f\n',sum(full(Vbi))); -fprintf('目标函数: %f\n',sum(full(PDbi))); +fprintf('目标函数: %f\n',sum(full(QDbi))); toc diff --git a/OPF_Init.m b/OPF_Init.m index f46528f..2ca898d 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -1,7 +1,7 @@ function [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) Loadi=find(QD~=0 | PD~=0); %Loadi=[1:Busnum]'; -RestraintCount=size(Loadi,1)*2+size(Loadi,1)+Busnum*2+Busnum+length(Loadi); %约束条件数,放开所有QD +RestraintCount=size(Loadi,1)*2+size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2; %约束条件数,放开所有QD t_Bal_volt=Volt(Balance); Volt=sparse(1*ones(1,Busnum)); Volt(Balance)=t_Bal_volt; diff --git a/func_ddh.m b/func_ddh.m index a4390b5..97823d0 100644 --- a/func_ddh.m +++ b/func_ddh.m @@ -51,7 +51,7 @@ t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT; sizeLoadi=size(Loadi,1)*2; ddh=[ sparse(sizeLoadi,ContrlCount); - sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum+length(Loadi)); - sparse(Busnum+length(Loadi),ContrlCount); + sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum+length(Loadi)*2); + sparse(Busnum+length(Loadi)*2,ContrlCount); ]; end \ No newline at end of file diff --git a/func_deltF.m b/func_deltF.m index 8933558..96ee1d3 100644 --- a/func_deltF.m +++ b/func_deltF.m @@ -6,5 +6,6 @@ deltF=[sparse(length(Loadi),1); sparse(2*Busnum,1); sparse(ones(Busnum,1)); sparse(ones(length(Loadi),1)); + sparse(ones(length(Loadi),1)); ]; end \ No newline at end of file diff --git a/func_deltG.m b/func_deltG.m index 2e7319c..631ea88 100644 --- a/func_deltG.m +++ b/func_deltG.m @@ -4,52 +4,75 @@ sizeLoadi=size(Loadi,1); dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1)); %dg32 PD+M*B+t+PD0 dg32_dPD=dg3_dPD; +%dg4 QD-M*B-t-QD0 dg4_dPD=sparse(size(Loadi,1),length(Loadi)); +%dg42 QD+M*B+t+QD0 +dg42_dPD=sparse(size(Loadi,1),length(Loadi)); dg5_dPD=sparse(size(Loadi,1),Busnum); dg6_dPD=dg5_dPD; %dg7 是 Vbi的约束 dg7_dPD=sparse(sizeLoadi,Busnum); dpdbi_dPD=sparse(sizeLoadi,sizeLoadi); +dqdbi_dPD=sparse(sizeLoadi,sizeLoadi); %% dg3_dQD=sparse(length(Loadi),length(Loadi)); dg32_dQD=sparse(sizeLoadi,sizeLoadi); dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1)); +dg42_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1)); dg5_dQD=sparse(size(Loadi,1),Busnum); dg6_dQD=dg5_dQD; dg7_dQD=sparse(sizeLoadi,Busnum); dpdbi_dQD=sparse(sizeLoadi,sizeLoadi); +dqdbi_dQD=sparse(sizeLoadi,sizeLoadi); %% dg3_dx=sparse(2*Busnum,sizeLoadi); dg32_dx=sparse(2*Busnum,sizeLoadi); dg4_dx=sparse(2*Busnum,length(Loadi)); +dg42_dx=sparse(2*Busnum,sizeLoadi); dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum); sparse(Busnum,Busnum); ]; dg6_dx=dg5_dx; dg7_dx=sparse(2*Busnum,Busnum); dpdbi_dx=sparse(2*Busnum,sizeLoadi); +dqdbi_dx=sparse(2*Busnum,sizeLoadi); %% dg3_dvbi=sparse(Busnum,sizeLoadi); dg32_dvbi=sparse(Busnum,sizeLoadi); dg4_dvbi=sparse(Busnum,length(Loadi)); +dg42_dvbi=sparse(Busnum,sizeLoadi); dg5_dvbi=sparse(-eye(Busnum,Busnum)); dg6_dvbi=sparse(eye(Busnum,Busnum)); dg7_dvbi=sparse(eye(Busnum,Busnum)); dpdbi_dvbi=sparse(Busnum,sizeLoadi); +dqdbi_dvbi=sparse(Busnum,sizeLoadi); %% dg3_dpdbi=sparse(-eye(sizeLoadi,sizeLoadi)); dg32_dpdbi=sparse(eye(sizeLoadi,sizeLoadi)); dg4_dpdbi=sparse(sizeLoadi,sizeLoadi); +dg42_dpdbi=sparse(sizeLoadi,sizeLoadi); dg5_dpdbi=sparse(sizeLoadi,Busnum); dg6_dpdbi=sparse(sizeLoadi,Busnum); dg7_dpdbi=sparse(sizeLoadi,Busnum); dpdbi_dpdbi=sparse(eye(sizeLoadi,sizeLoadi)); +dqdbi_dpdbi=sparse(sizeLoadi,sizeLoadi); %% -deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg5_dPD,dg6_dPD,dg7_dPD,dpdbi_dPD; - dg3_dQD,dg32_dQD,dg4_dQD,dg5_dQD,dg6_dQD,dg7_dQD,dpdbi_dQD; - dg3_dx,dg32_dx,dg4_dx,dg5_dx,dg6_dx,dg7_dx,dpdbi_dx; - dg3_dvbi,dg32_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi,dpdbi_dvbi; - dg3_dpdbi,dg32_dpdbi,dg4_dpdbi,dg5_dpdbi,dg6_dpdbi,dg7_dpdbi,dpdbi_dpdbi; +dg3_dqdbi=sparse(sizeLoadi,sizeLoadi); +dg32_dqdbi=sparse(sizeLoadi,sizeLoadi); +dg4_dqdbi=sparse(-eye(sizeLoadi,sizeLoadi)); +dg42_dqdbi=sparse(eye(sizeLoadi,sizeLoadi)); +dg5_dqdbi=sparse(sizeLoadi,Busnum); +dg6_dqdbi=sparse(sizeLoadi,Busnum); +dg7_dqdbi=sparse(sizeLoadi,Busnum); +dpdbi_dqdbi=sparse(sizeLoadi,sizeLoadi); +dqdbi_dqdbi=sparse(eye(sizeLoadi,sizeLoadi)); +%% +deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg42_dPD,dg5_dPD,dg6_dPD,dg7_dPD,dpdbi_dPD,dqdbi_dPD; + dg3_dQD,dg32_dQD,dg4_dQD,dg42_dQD,dg5_dQD,dg6_dQD,dg7_dQD,dpdbi_dQD,dqdbi_dQD; + dg3_dx,dg32_dx,dg4_dx,dg42_dx,dg5_dx,dg6_dx,dg7_dx,dpdbi_dx,dqdbi_dx; + dg3_dvbi,dg32_dvbi,dg4_dvbi,dg42_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi,dpdbi_dvbi,dqdbi_dvbi; + dg3_dpdbi,dg32_dpdbi,dg4_dpdbi,dg42_dpdbi,dg5_dpdbi,dg6_dpdbi,dg7_dpdbi,dpdbi_dpdbi,dqdbi_dpdbi; + dg3_dqdbi,dg32_dqdbi,dg4_dqdbi,dg42_dqdbi,dg5_dqdbi,dg6_dqdbi,dg7_dqdbi,dpdbi_dqdbi,dqdbi_dqdbi; ]; % deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg5_dPD,dg6_dPD,dg7_dPD; % dg3_dQD,dg32_dQD,dg4_dQD,dg5_dQD,dg6_dQD,dg7_dQD; diff --git a/func_deltH.m b/func_deltH.m index cddcb8e..ae4443d 100644 --- a/func_deltH.m +++ b/func_deltH.m @@ -4,6 +4,7 @@ dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loa dH_dx = jacobian_M(Busnum,Volt,Y,Angle,UAngel,r,c); %形成雅克比矩阵 dH_dvbi=sparse(Busnum,2*Busnum); dH_dpdbi=sparse(length(Loadi),2*Busnum); -deltH=[dH_dPD;dH_dQD;dH_dx';dH_dvbi;dH_dpdbi]; +dH_dqdbi=sparse(length(Loadi),2*Busnum); +deltH=[dH_dPD;dH_dQD;dH_dx';dH_dvbi;dH_dpdbi;dH_dqdbi]; end \ No newline at end of file