From f9dc83b8626e09e659e8ee49d129a2ea70b4484d Mon Sep 17 00:00:00 2001 From: facat Date: Wed, 6 Jun 2012 10:19:08 +0800 Subject: [PATCH] =?UTF-8?q?=E9=9A=8F=E6=9C=BA=E6=B5=8B=E8=AF=95=E7=BC=BA?= =?UTF-8?q?=E8=B4=9F=E8=8D=B7?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- FormH.m | 6 ++-- FormLw.asv | 12 ++++--- FormLw.m | 12 ++++--- FormLz.m | 12 ++++--- FormYY.m | 2 +- OPF.asv | 22 ++++++++---- OPF.m | 23 ++++++++++--- OPF_Init.m | 12 +++---- SolveIt.m | 5 +-- func_ddh.m | 9 +++-- func_deltF.m | 2 +- test.m | 97 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 12 files changed, 172 insertions(+), 42 deletions(-) create mode 100644 test.m diff --git a/FormH.m b/FormH.m index 5b735eb..6d5cb72 100644 --- a/FormH.m +++ b/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'); -% QD(QD~=0)=PD(QD~=0)./tan(QDcos); -% QD(QD_NON_ZERO_IND)=QD_NON_ZERO; +QDcos=textread('D:\Project\最小化潮流\最小潮流算例\仙海919PDQDglys.txt'); +QD(QD~=0)=PD(QD~=0)./tan(QDcos); +%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'; diff --git a/FormLw.asv b/FormLw.asv index 1c9e4b4..5dd6cf0 100644 --- a/FormLw.asv +++ b/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)=0.8*PDU(PDU<0); -PDU(PDU==0)=.2; +% PDU(PDU>0)=1200*PDU(PDU>0); +% PDU(PDU<0)=-800*PDU(PDU<0); +% PDU(PDU==0)=200; +PDU=1000000*ones(length(Loadi),1); t1=([PU',QU',PDU',VoltU])'; t2=Mat_G+Init_U'-t1; Lw=t2; diff --git a/FormLw.m b/FormLw.m index 37351ca..e7ecda5 100644 --- a/FormLw.m +++ b/FormLw.m @@ -1,12 +1,14 @@ function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK) -%KK=999; -PU=1*GenU(:,2);%发电机有功上界 +KK=999; +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)=.80*PDU(PDU<0); -PDU(PDU==0)=.20; +PDU(PDU>0)=1.200*PDU(PDU>0); +PDU(PDU<0)=0.800*PDU(PDU<0); +PDU(PDU==0)=0.400; +%PDU=10*ones(length(Loadi),1); t1=([PU',QU',PDU',VoltU])'; t2=Mat_G+Init_U'-t1; Lw=t2; diff --git a/FormLz.m b/FormLz.m index eedd55c..1964285 100644 --- a/FormLz.m +++ b/FormLz.m @@ -1,12 +1,14 @@ function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK) -%KK=999; -PL=1*GenL(:,2);%发电机有功下界 +KK=999; +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)=1.20*PDL(PDL<0); -PDL(PDL==0)=-.20; +PDL(PDL>0)=0.800*PDL(PDL>0); +PDL(PDL<0)=1.200*PDL(PDL<0); +PDL(PDL==0)=-0.400; +%PDL=-10*ones(length(Loadi),1); t1=([PL',QL',PDL',VoltL])'; t2=Mat_G-Init_L'-t1; Lz=t2; diff --git a/FormYY.m b/FormYY.m index 0049526..cbad106 100644 --- a/FormYY.m +++ b/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; diff --git a/OPF.asv b/OPF.asv index a5d5ebb..bc967f6 100644 --- a/OPF.asv +++ b/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); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,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,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'); -% PD=full(PD); +Volt=full(Volt'); +PD=full(PD); +%% 统计PD误差 +abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) ) toc diff --git a/OPF.m b/OPF.m index 18613dc..c5802c7 100644 --- a/OPF.m +++ b/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'); -%PVi电压节点序号 -%PVu电压节点电压标幺值 +[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'); + +%% 计算功率因数 +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'); -% PD=full(PD); +%Volt=full(Volt'); +%PD=full(PD); +%% 统计PD误差 +maxPDError=max(abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) )); toc diff --git a/OPF_Init.m b/OPF_Init.m index 1f04e02..0ca0a4f 100644 --- a/OPF_Init.m +++ b/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); - -wD=ones(size(Loadi,1),1); +randPDind=randInt(1:HH); +wD=0.1*ones(size(Loadi,1),1); wD(randPDind)=0;%一些负荷不约束 -%wD(12)=0; +%wD(7)=0; +% wD(11)=0; PD=1*PD0; end \ No newline at end of file diff --git a/SolveIt.m b/SolveIt.m index 672f754..23d81cc 100644 --- a/SolveIt.m +++ b/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'); -t1=diag(Init_L.\Init_Z-Init_U.\Init_W); +H=-deltdeltF+ddh; +%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; diff --git a/func_ddh.m b/func_ddh.m index 22ed8fb..e678761 100644 --- a/func_ddh.m +++ b/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 ; - ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT; +% t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV; +% ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT; +% ]; +t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT; + ddPdTdV+ddQdTdV,ddPdTdT+ddQdTdT; ]; sizePGi=size(PGi,1); sizePVi=size(PVi,1); diff --git a/func_deltF.m b/func_deltF.m index 3cf4c06..3217b2f 100644 --- a/func_deltF.m +++ b/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); ]; diff --git a/test.m b/test.m new file mode 100644 index 0000000..e7b3690 --- /dev/null +++ b/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 + +