From 17412d8ec53f0a6ec01307d303a779087e7aafea Mon Sep 17 00:00:00 2001 From: "dugg@lab-desk" Date: Sun, 19 Oct 2014 11:42:44 +0800 Subject: [PATCH] =?UTF-8?q?=E6=8A=8A=E5=8E=9F=E6=9D=A5=E7=9A=84OPF?= =?UTF-8?q?=E9=87=8D=E6=96=B0=E5=8A=A0=E5=85=A5=E8=BF=99=E9=87=8C=E9=9D=A2?= =?UTF-8?q?=EF=BC=8C=E6=94=B6=E6=95=9B=E4=BA=86=EF=BC=8C=E4=BD=86=E6=80=BB?= =?UTF-8?q?=E8=A7=89=E5=BE=97=E6=9C=89=E5=9C=B0=E6=96=B9=E4=B8=8D=E5=A4=AA?= =?UTF-8?q?=E5=AF=B9=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dugg@lab-desk --- FormG.m | 14 +++++++++---- FormH.m | 14 +++---------- FormLw.m | 9 +-------- FormLz.m | 11 +---------- IPMLoop.m | 49 +++++++++++++++++++++++++--------------------- Modification.m | 27 ++++++++++++++----------- SolveIt.m | 25 ++++++++++------------- feeder13/data1.txt | 2 +- func_ddg.m | 2 +- func_deltF.m | 18 ++++++++--------- func_deltG.m | 15 -------------- func_deltH.m | 4 ---- func_deltdeltF.m | 15 -------------- 13 files changed, 79 insertions(+), 126 deletions(-) diff --git a/FormG.m b/FormG.m index 3325803..763ef96 100644 --- a/FormG.m +++ b/FormG.m @@ -1,8 +1,14 @@ function Mat_G=FormG(I1r,I1i) -% Mat_G=[ -% I1r; -% I1i; -% ]; +%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=[ I1r.^2+I1i.^2; ]; diff --git a/FormH.m b/FormH.m index 9aff295..4431224 100644 --- a/FormH.m +++ b/FormH.m @@ -1,15 +1,6 @@ function Mat_H=FormH(fsY11,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance) -%% -%QDcos=textread('D:\Project\最小化潮流\最小潮流算例\仙海919PDQDglys.txt'); -%QD(QD~=0)=PD(QD~=0)./tan(QDcos); -%QD(QD_NON_ZERO_IND)=QD_NON_ZERO; -%% -%PD(Loadi)=QD(Loadi)./tan(acos(0.98)); -% AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); -% dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt; -% dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt; -% -% Mat_H=[dP;dQ;]; + +%%%%一下是学姐给的公式 H=[ real(fsY11),-imag(fsY11); imag(fsY11),real(fsY11); @@ -21,4 +12,5 @@ Mat_H=H*[V1r;V1i]; Mat_H=Mat_H+sparse(Loadi,1,-I1r,busNum*2,1)+sparse(Loadi+busNum,1,-I1i,busNum*2,1); Mat_H(Balance)=Mat_H(Balance)-BalI1r; Mat_H(Balance+busNum)=Mat_H(Balance+busNum)-BalI1i; + end \ No newline at end of file diff --git a/FormLw.m b/FormLw.m index 168d52b..c036f78 100644 --- a/FormLw.m +++ b/FormLw.m @@ -1,12 +1,5 @@ function Lw=FormLw(Loadi,Mat_G,Init_U) -% VoltU=sparse(1.2*ones(Busnum,1)); -% PDU=sparse(2*ones(length(Loadi),1)); -% QDU=sparse(2*ones(length(Loadi),1)); -% t1=[PDU; -% QDU; -% VoltU]; -% t2=Mat_G+Init_U-t1; -% Lw=t2; + upper=2*sparse(ones(length(Loadi)*1,1)); Lw=Mat_G+Init_U-upper; end \ No newline at end of file diff --git a/FormLz.m b/FormLz.m index 05c019e..992d38f 100644 --- a/FormLz.m +++ b/FormLz.m @@ -1,14 +1,5 @@ function Lz=FormLz(Loadi,Mat_G,Init_L) -% VoltL=sparse(0.8*ones(Busnum,1)); -% PDL=sparse(-2*ones(length(Loadi),1)); -% QDL=sparse(-2*ones(length(Loadi),1)); -% t1=[PDL; -% QDL; -% VoltL -% ]; -% t2=Mat_G-Init_L-t1; -% Lz=t2; + lower=0*sparse(ones(length(Loadi)*1,1)); Lz=Mat_G-Init_L-lower; -% Lz=-[] end \ No newline at end of file diff --git a/IPMLoop.m b/IPMLoop.m index 3a4e3e3..aed3020 100644 --- a/IPMLoop.m +++ b/IPMLoop.m @@ -2,8 +2,8 @@ function [ V1r,V1i,I1r,I1i ] = IPMLoop(measurement,BalI1r,BalI1i,busNum,Loadi,fs %把每个序的循环写在这个函数中。其实也就是内点法循环。 V1r=1*ones(busNum,1); V1i=0*ones(busNum,1); -I1r=0.01*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反 -I1i=0.01*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反 +I1r=0.1*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反 +I1i=0.1*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反 KK=0; plotGap=zeros(1,60); %初始化 @@ -12,32 +12,36 @@ RestraintCount=length(Loadi)*1; ContrlCount=busNum*2+length(Loadi)*2; CenterA=0.1; Init_Z=sparse(ones(RestraintCount,1)); -Init_W=sparse(-1*ones(RestraintCount,1)); +Init_W=sparse(-1.0*ones(RestraintCount,1)); Init_L=1*sparse(ones(RestraintCount,1)); Init_U=1*sparse(ones(RestraintCount,1)); Init_Y=sparse(2*busNum,1);%等式约束乘子 Gap=(Init_L'*Init_Z-Init_U'*Init_W); -% PG=sparse(busNum,1); -% PG(Balance)=0.1105; -% QG=sparse(busNum,1); -% QG(Balance)=0.0984; -% SEPD(2)=0.1105; -% SEQD(2)=0.0984; -while Gap>1e-5 && KK<20 - KK=KK+1; +kmax=60; +while(abs(Gap)>0.000001) + if KK>kmax + break; + end + plotGap(KK+1)=Gap; Init_u=Gap/2/RestraintCount*CenterA; + %% 开始计算OPF + %% 形成等式约束的雅克比 deltH=func_deltH(busNum,fsY1,Loadi,Balance); + %% 形成不等式约束的雅克比 deltG=func_deltG(busNum,Loadi,I1r,I1i); -% deltG=0; -% L_1Z=diag(Init_Z./Init_L); -% U_1W=diag(Init_W./Init_U); + %% + L_1Z=diag(Init_Z./Init_L); + U_1W=diag(Init_W./Init_U); + %% 形成海森阵 deltdeltF=func_deltdeltF(busNum,fsY1,Loadi); -% deltdeltF=0; -% ddh=func_ddh(SEVmf1,Init_Y,busNum,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi,ContrlCount); + %% 形成ddHy +% ddh=func_ddh(busNum,Loadi,Init_Z,Init_W); ddh=0; + %% 开始构建ddg ddg=func_ddg(busNum,Loadi,Init_Z,Init_W); + %% 开始构建deltF deltF=func_deltF(measurement,busNum,fsY1,Loadi,V1r,V1i,I1r,I1i); -% deltF=0; + %% Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1); Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1); Mat_G=FormG(I1r,I1i); @@ -46,14 +50,15 @@ while Gap>1e-5 && KK<20 Lz=FormLz(Loadi,Mat_G,Init_L); Lw=FormLw(Loadi,Mat_G,Init_U); Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); -% YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); + YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); %% 开始解方程 - plotGap(KK)=Gap; - fprintf('迭代次数 %d Gap %f\n',KK,plotGap(KK)); - XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,busNum,Loadi); + 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,busNum); + %%取各分量 [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum); - [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi,Vref); + [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi); Gap=(Init_L'*Init_Z-Init_U'*Init_W); + fprintf('Gap %f\n',full(Gap)); + KK=KK+1; end f=sum(([real(measurement);imag(measurement)]-[-I1r;-I1i]).^2); diff --git a/Modification.m b/Modification.m index fd47839..a1a94e5 100644 --- a/Modification.m +++ b/Modification.m @@ -1,22 +1,27 @@ -function [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY, ... - V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi,Vref) +function [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,Busnum,Loadi) AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU); +% fprintf('AlphaP %f\n',full(AlphaP)); AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW); +% fprintf('AlphaD %f\n',full(AlphaD)); + Init_Z=Init_Z+AlphaD*deltZ; Init_L=Init_L+AlphaP*deltL; Init_W=Init_W+AlphaD*deltW; Init_U=Init_U+AlphaP*deltU; Init_Y=Init_Y+AlphaD*deltY; -t=deltX(1:busNum*2); -V1r=V1r+AlphaP*t(1:busNum); -V1r(Balance)=Vref; -V1i=V1i+AlphaP*t(busNum+1:end); +%PG(PGi)=PG(PGi)+deltX(1:size(PGi,1)); +% PG(PGi)=PG(PGi)+AlphaP*deltX(1:size(PGi,1)); +V1r=V1r+AlphaP*deltX(1:Busnum); +V1r(Balance)=1; +V1i=V1i+AlphaP*deltX(Busnum+1:2*Busnum); V1i(Balance)=0; -t=deltX(busNum*2+1:ContrlCount); +%QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) ); +% QG(PVi)=QG(PVi)+AlphaP*deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) ); +t=deltX(Busnum*2+1:ContrlCount); I1r=I1r+AlphaP*t(1:length(Loadi)); I1i=I1i+AlphaP*t(length(Loadi)+1:end); -%balVolt=1;%Volt(Balance); -%Volt=Volt+AlphaP*t(1:Busnum); -%Volt(Balance)=balVolt; -%UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum); + +%Volt=Volt+AlphaP*t(2:2:2*Busnum);暂时改一下20111227 +%UAngel=UAngel+AlphaP*t(1:2:2*Busnum);暂时改一下20111227 + end \ No newline at end of file diff --git a/SolveIt.m b/SolveIt.m index a7d573a..f1e21e4 100644 --- a/SolveIt.m +++ b/SolveIt.m @@ -1,31 +1,26 @@ -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,Lx,Balance,Busnum,Loadi) +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,ConstraintCount,Lx,Balance,Busnum) LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx); -H=-deltdeltF+ddh+ddg; -t1=diag(Init_Z./Init_L-Init_W./Init_U); +H=-deltdeltF+ddh+ddg;%+ddg*(Init_Z'+Init_W'); +t1=diag(Init_L.\Init_Z-Init_U.\Init_W); t2=-deltG*( t1 )*deltG'; aa=[ (H+t2),deltH; - deltH',zeros(length(Init_Y)); + deltH',sparse(length(Init_Y),length(Init_Y)); ]; yy=[LxComa;-Ly]; -%% 平衡节点电压实部 +% t=size(PVi,1)+size(PGi,1);%电压不变 +% yy(t+Balance)=0; aa(Balance,:)=0; aa(:,Balance)=0; aa=aa+sparse(Balance,Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); -deltG(Balance,:)=0; -%% 平衡节点电压虚部 -t=Busnum*1; + +%暂时改一下 +t=Busnum;%相角不变 aa(t+Balance,:)=0; aa(:,t+Balance)=0; -aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); -deltG(t+Balance,:)=0; -%% +aa(t+Balance,t+Balance)=1; dxdy=aa\yy; -%% KLU -%spy(aa) -%dxdy = klu(aa,'\',full(yy)); -%% dX=dxdy(1:ContrlCount); dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum); dL=Lz+deltG'*dX; diff --git a/feeder13/data1.txt b/feeder13/data1.txt index 0b1c5b9..33f02f4 100644 --- a/feeder13/data1.txt +++ b/feeder13/data1.txt @@ -5,8 +5,8 @@ 0 632 30 0 0 0 -632 120 120 120 110 90 90 631 90 100 110 80 90 90 +632 120 120 120 110 90 90 0 diff --git a/func_ddg.m b/func_ddg.m index f79563d..68ba20f 100644 --- a/func_ddg.m +++ b/func_ddg.m @@ -1,5 +1,5 @@ function ddg=func_ddg(busNum,Loadi,Init_Z,Init_W) -% ddg=0; + Init_ZW=Init_Z+Init_W; t=[2*diag( ones(length(Loadi),1).*Init_ZW),zeros(length(Loadi)); zeros(length(Loadi)),2*diag( ones(length(Loadi),1).*Init_ZW); diff --git a/func_deltF.m b/func_deltF.m index 980e4c7..3594b25 100644 --- a/func_deltF.m +++ b/func_deltF.m @@ -1,15 +1,15 @@ function deltF=func_deltF(measurement,busNum,fsY11,Loadi,V1r,V1i,I1r,I1i) -% H=[ -% real(fsY11),-imag(fsY11),zeros(busNum,length(Loadi)*2); -% imag(fsY11),real(fsY11),zeros(busNum,length(Loadi)*2); -% zeros(length(Loadi),busNum*2), eye(length(Loadi)),zeros(length(Loadi)); -% zeros(length(Loadi),busNum*2), zeros(length(Loadi)),eye(length(Loadi)); -% ]; -% X=[V1r;V1i;I1r;I1i]; -% deltF=-H'*(measurement-H*X); +%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])';%暂时用土办法处理一下 +%% deltF=[ zeros(busNum*2,1); -2*( [real(measurement);imag(measurement)]-[I1r;I1i]); ]; -% deltF=0; + end \ No newline at end of file diff --git a/func_deltG.m b/func_deltG.m index 0a8756f..97bee74 100644 --- a/func_deltG.m +++ b/func_deltG.m @@ -1,19 +1,4 @@ function deltG=func_deltG(busNum,Loadi,I1r,I1i) -% dgV_dPD=sparse(length(Loadi),busNum); -% dgV_dQD=sparse(length(Loadi),busNum); -% dgV_dx=[sparse(1:busNum,1:busNum,ones(busNum,1),busNum,busNum); -% sparse(busNum,busNum); -% ]; -% % -% dgPD_dPD=sparse(eye(length(Loadi))); -% dgPD_dQD=sparse(length(Loadi),length(Loadi)); -% dgPD_dx=sparse(busNum*2,length(Loadi)); -% % -% dgQD_dPD=sparse(length(Loadi),length(Loadi)); -% dgQD_dQD=sparse(eye(length(Loadi))); -% dgQD_dx=sparse(busNum*2,length(Loadi)); -% %% - deltG=[ zeros(busNum*2,length(Loadi)*1); 2*eye(length(Loadi)*1)*diag(I1r); diff --git a/func_deltH.m b/func_deltH.m index c69b752..bff54ef 100644 --- a/func_deltH.m +++ b/func_deltH.m @@ -1,8 +1,4 @@ function deltH=func_deltH(busNum,fsY11,Loadi,Balance) -% dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)]; -% dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)]; -% dH_dx = jacobian(Busnum,Volt,Y,Angle,UAngel,r,c); %形成雅克比矩阵 -% deltH=[dH_dPD;dH_dQD;dH_dx']; H=[ real(fsY11),-imag(fsY11); imag(fsY11),real(fsY11); diff --git a/func_deltdeltF.m b/func_deltdeltF.m index 2003458..d07d26b 100644 --- a/func_deltdeltF.m +++ b/func_deltdeltF.m @@ -1,22 +1,7 @@ function deltdeltF=func_deltdeltF(busNum,fsY11,Loadi) -% H=[ -% real(fsY11),-imag(fsY11),zeros(busNum,length(Loadi)*2); -% imag(fsY11),real(fsY11),zeros(busNum,length(Loadi)*2); -% zeros(length(Loadi),busNum*2), eye(length(Loadi)),zeros(length(Loadi)); -% zeros(length(Loadi),busNum*2), zeros(length(Loadi)),eye(length(Loadi)); -% ]; -% H=[ -% real(fsY11),-imag(fsY11),zeros(busNum,length(Loadi)*2); -% imag(fsY11),real(fsY11),zeros(busNum,length(Loadi)*2); -% zeros(length(Loadi),busNum*2), eye(length(Loadi)),zeros(length(Loadi)); -% zeros(length(Loadi),busNum*2), zeros(length(Loadi)),eye(length(Loadi)); -% ]; - -% deltdeltF=H'*H; deltdeltF=[ zeros(busNum*2,busNum*2+length(Loadi)*2); zeros(length(Loadi)*2,busNum*2),2*eye(length(Loadi)*2); ]; -% deltdeltF=0; end \ No newline at end of file