facat
/
pes2014-debug2-laplace
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

加入部分负荷未知的情况。修复了有功无功上下界的一个bug。 

Signed-off-by: facat <dmy@dmy-PC.(none)>
This commit is contained in:
facat 2012-05-25 15:38:31 +08:00
parent f2781393a7
commit 460c49829b
19 changed files with 122 additions and 136 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
FormG.m
View File

@ -1,4 +1,4 @@
function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD) function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,Loadi)
%t1=PG(PVi); %t1=PG(PVi);
%GP=t1;%发电机P %GP=t1;%发电机P
%GP=[4.5 4.5]'; %GP=[4.5 4.5]';
@ -13,7 +13,7 @@ Mat_G=[
%GP; %GP;
PG(PGi); PG(PGi);
QG(PVi); QG(PVi);
sparse(PD); sparse(PD(Loadi));
%GQ(PVi); %GQ(PVi);
%[0 1.45]'; %[0 1.45]';
Volt'; Volt';

1
FormH.m
View File

@ -5,7 +5,6 @@ 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); AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
%dP=PG-PD-diag(Volt)*Y*cos(AngleIJ)*Volt';
dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt'; dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt'; dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';

10
FormLw.m
View File

@ -1,14 +1,14 @@
function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0) function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi)
PU=1*GenU(:,2);%发电机有功上界 PU=1*GenU(:,2);%发电机有功上界
QU=PVQU(:,1);%发电机无功上界 QU=PVQU(:,1);%发电机无功上界
VoltU=1.1*ones(1,Busnum); VoltU=1.1*ones(1,Busnum);
PDU=PD0; PDU=PD0(Loadi);
PDU(PD0>0)=1.2*PD0(PD0>0); PDU(PDU>0)=1.2*PDU(PDU>0);
PDU(PDU<0)=0.8*PDU(PDU<0);
PDU(PDU==0)=.2; PDU(PDU==0)=.2;
PDU(PDU<0)=0.8*PD0(PD0<0);
%PDU=25*ones(Busnum,1)';
t1=([PU',QU',PDU',VoltU])'; t1=([PU',QU',PDU',VoltU])';
t2=Mat_G+Init_U'-t1; t2=Mat_G+Init_U'-t1;
Lw=t2; Lw=t2;
end end

24
FormLz.asv
View File

@ -1,10 +1,28 @@
function Lz=FormLz(Mat_G,Init_L,GenL,LinePLimt,LineCount,Busnum,PVQL) function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi)
PL=GenL(:,2);%发电机有功下界 PL=1*GenL(:,2);%发电机有功下界
QL=PVQL(:,1);%发电机无功下界 QL=PVQL(:,1);%发电机无功下界
VoltL=0.9*ones(1,Busnum); VoltL=0.9*ones(1,Busnum);
t1=([PL',QL',VoltL])'; PDL=PD0(Loadi);
PDL(PD0>0)=0.8*PD0(PDL>0);
PDL(PD0==0)=-.2;
PDL(PDL<0)=1.2*PD0(PDL<0);
t1=([PL',QL',PDL',VoltL])';
t2=Mat_G-Init_L'-t1; t2=Mat_G-Init_L'-t1;
Lz=t2; Lz=t2;
% PL=1*GenL(:,2);%发电机有功下界
% QL=PVQL(:,1);%发电机无功下界
% VoltL=0.9*ones(1,Busnum);
% PDL=PD0;
% PDL(PD0>0)=0.8*PD0(PD0>0);
% PDL(PD0==0)=-.2;
% PDL(PD0<0)=1.2*PD0(PD0<0);
%
% t1=([PL',QL',PDL',VoltL])';
%
% t2=Mat_G-Init_L'-t1;
% Lz=t2;
end end

13
FormLz.m
View File

@ -1,15 +1,14 @@
function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0) function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi)
PL=1*GenL(:,2);%发电机有功下界 PL=1*GenL(:,2);%发电机有功下界
QL=PVQL(:,1);%发电机无功下界 QL=PVQL(:,1);%发电机无功下界
VoltL=0.9*ones(1,Busnum); VoltL=0.9*ones(1,Busnum);
PDL=PD0; PDL=PD0(Loadi);
PDL(PD0>0)=0.8*PD0(PD0>0); PDL(PDL>0)=0.8*PDL(PDL>0);
PDL(PD0==0)=-.2; PDL(PDL<0)=1.2*PDL(PDL<0);
PDL(PD0<0)=1.2*PD0(PD0<0); PDL(PDL==0)=-.2;
t1=([PL',QL',PDL',VoltL])'; t1=([PL',QL',PDL',VoltL])';
t2=Mat_G-Init_L'-t1; t2=Mat_G-Init_L'-t1;
Lz=t2; Lz=t2;
end end

8
Modification.m
View File

@ -1,4 +1,4 @@
function [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) function [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)
AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU); AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
fprintf('AlphaP %f\n',full(AlphaP)); fprintf('AlphaP %f\n',full(AlphaP));
AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW); AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
@ -13,9 +13,9 @@ Init_Y=Init_Y+AlphaD*deltY';
PG(PGi)=PG(PGi)+AlphaP*deltX(1:size(PGi,1)); PG(PGi)=PG(PGi)+AlphaP*deltX(1:size(PGi,1));
%QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) ); %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) ); QG(PVi)=QG(PVi)+AlphaP*deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
t=deltX(size(PVi,1)+size(PGi,1)+1:size(PVi,1)+size(PGi,1)+Busnum); t=deltX(size(PVi,1)+size(PGi,1)+1:size(PVi,1)+size(PGi,1)+size(Loadi,1));
PD=PD+AlphaP*t; PD(Loadi)=PD(Loadi)+AlphaP*t;
t=deltX(size(PVi,1)+size(PGi,1)+Busnum+1:ContrlCount)'; t=deltX(size(PVi,1)+size(PGi,1)+size(Loadi,1)+1:ContrlCount)';
t(Busnum+Balance)=0; t(Busnum+Balance)=0;
%Volt=Volt+AlphaP*t(2:2:2*Busnum);暂时改一下20111227 %Volt=Volt+AlphaP*t(2:2:2*Busnum);暂时改一下20111227
%UAngel=UAngel+AlphaP*t(1:2:2*Busnum);暂时改一下20111227 %UAngel=UAngel+AlphaP*t(1:2:2*Busnum);暂时改一下20111227

15
OPF.asv
View File

@ -1,18 +1,14 @@
tic tic
clear clear
%[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,LineLimti,LineLimtj,LinePLimt,PG,QG,PD,QD,CenterA,LineCount,PGi,PVQU,PVQL]=pf('5sj.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('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('ieee3001PG.dat');
GB=full(GB);
%PVi电压节点序号 %PVi电压节点序号
%PVu电压节点电压标幺值 %PVu电压节点电压标幺值
Volt; Volt;
UAngel*180/3.1415926; UAngel*180/3.1415926;
%sprintf('%f\n',Volt);
%sprintf('%f\n',Angel);
%% 初值 %% 初值
PG0=PG; PG0=PG;
PD0=PD; PD0=PD;
[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD]=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,randPDind]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0);
Gap=(Init_L*Init_Z'-Init_U*Init_W'); Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0; KK=0;
plotGap=zeros(1,50); plotGap=zeros(1,50);
@ -28,10 +24,10 @@ while(abs(Gap)>Precision)
end end
plotGap(KK+1)=Gap; plotGap(KK+1)=Gap;
Init_u=Gap/2/RestraintCount*CenterA; Init_u=Gap/2/RestraintCount*CenterA;
AngleIJMat=repmat(UAngel',1,Busnum)-repmat(UAngel,Busnum,1); AngleIJMat=0;
%% 开始计算OPF %% 开始计算OPF
%% 形成等式约束的雅克比 %% 形成等式约束的雅克比
deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi); deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle);
%% 形成不等式约束的雅克比 %% 形成不等式约束的雅克比
deltG=func_deltG(Busnum,PVi,PGi); deltG=func_deltG(Busnum,PVi,PGi);
%% %%
@ -40,7 +36,7 @@ while(abs(Gap)>Precision)
%% 形成海森阵 %% 形成海森阵
deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0); deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0);
%% 形成ddHy %% 形成ddHy
ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y); ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle);
%% 开始构建ddg %% 开始构建ddg
ddg=func_ddg(PGi,PVi,Busnum,RestraintCount); ddg=func_ddg(PGi,PVi,Busnum,RestraintCount);
%% 开始构建deltF %% 开始构建deltF
@ -50,7 +46,7 @@ while(abs(Gap)>Precision)
Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1); Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
Lul=Init_L'.*Init_Z'-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);
Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,QD_NON_ZERO,QD_NON_ZERO_IND); 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; Ly=Mat_H;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0); Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0); Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0);
@ -65,7 +61,6 @@ while(abs(Gap)>Precision)
KK=KK+1; KK=KK+1;
end end
fprintf('迭代次数%d\n',KK); fprintf('迭代次数%d\n',KK);
%CalCost(GenC,PG,PGi);
ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD) ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD)
DrawGap(plotGap); DrawGap(plotGap);
Volt=Volt'; Volt=Volt';

30
OPF.m
View File

@ -8,11 +8,11 @@ UAngel*180/3.1415926;
%% 初值 %% 初值
PG0=PG; PG0=PG;
PD0=PD; PD0=PD;
[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD]=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'); Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0; KK=0;
plotGap=zeros(1,50); 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; kmax=600;
%% 20120523 临时 %% 20120523 临时
QD_NON_ZERO=QD(PD==0 & QD~=0); QD_NON_ZERO=QD(PD==0 & QD~=0);
@ -27,42 +27,44 @@ while(abs(Gap)>Precision)
AngleIJMat=0; AngleIJMat=0;
%% 开始计算OPF %% 开始计算OPF
%% 形成等式约束的雅克比 %% 形成等式约束的雅克比
deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle); deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,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); L_1Z=diag(Init_Z./Init_L);
U_1W=diag(Init_W./Init_U); U_1W=diag(Init_W./Init_U);
%% 形成海森阵 %% 形成海森阵
deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0); deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0);
%% 形成ddHy %% 形成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
ddg=func_ddg(PGi,PVi,Busnum,RestraintCount); ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi);
%% 开始构建deltF %% 开始构建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); Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
Lul=Init_L'.*Init_Z'-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); 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; Ly=Mat_H;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0); Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0); Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx); 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); [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'); Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=KK+1; KK=KK+1;
end end
fprintf('迭代次数%d\n',KK); fprintf('迭代次数%d\n',KK);
ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD) ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
DrawGap(plotGap); DrawGap(plotGap);
Volt=Volt'; %%
Volt=full(Volt');
PD=full(PD);
toc toc

42
OPF_Init.asv
View File

@ -1,27 +1,29 @@
function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0) function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0)
RestraintCount=size(PVi,1)+size(PGi,1)+Busnum*2; %约束条件数 RestraintCount=size(PVi,1)+size(PGi,1)+Busnum*2; %约束条件数
t_Bal_volt=Volt(Balance); t_Bal_volt=Volt(Balance);
Volt=ones(1,Busnum); Volt=sparse(ones(1,Busnum));
Volt(Balance)=t_Bal_volt; Volt(Balance)=t_Bal_volt;
%Volt(Balance)=1; UAngel=sparse(1,Busnum);
UAngel=zeros(1,Busnum); Init_Z=sparse(ones(1,RestraintCount));
Init_Z=ones(1,RestraintCount); Init_W=sparse(-1*ones(1,RestraintCount));
Init_W=-1*ones(1,RestraintCount); Init_L=sparse(ones(1,RestraintCount));
Init_L=ones(1,RestraintCount); Init_U=sparse(ones(1,RestraintCount));
Init_U=ones(1,RestraintCount); Init_Y=sparse(1,2*Busnum);%与学姐一致
%Init_Y=zeros(1,2*Busnum); tPU=sparse(GenU(:,2));% 发电机有功上限
%Init_Y=ones(1,2*Busnum); tQU=sparse(PVQU(:,1));% 无功上限
Init_Y=zeros(1,2*Busnum);%与学姐一致 tPL=sparse(GenL(:,2));% 发电机有功下限
%Init_Y(1:2:2*Busnum)=1e-10; tQL=sparse(PVQL(:,1));% 无功下限
%Init_Y(2:2:2*Busnum)=-1e-10;
tPU=GenU(:,2);% 发电机有功上限
tQU=PVQU(:,1);% 无功上限
tPL=GenL(:,2);% 发电机有功下限
tQL=PVQL(:,1);% 无功下限
%PG(4:5)=[4.5 4.5];
PG(PGi)=(tPU+tPL)/2; PG(PGi)=(tPU+tPL)/2;
%QG(4:5)=[0 1.45];
QG(PVi)=(tQU+tQL)/2; QG(PVi)=(tQU+tQL)/2;
wG=ones(size(PGi,1),1);
randInt=randperm(300);
asc_randInt=sort(randInt,2,'ascend' );
randPDind=randInt(1:200);
wD=ones(Busnum,1); wD=ones(Busnum,1);
PD=.5*PD0 wD(randPDind)=0;%一些负荷不约束
%wD(Balance)=0;
PD0(randPDind)=0;
PD=1*PD0;
%PD(PD==0)=.2;
end end

21
OPF_Init.m
View File

@ -1,28 +1,33 @@
function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0) 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)
RestraintCount=size(PVi,1)+size(PGi,1)+Busnum*2; %Ô¼ÊøÌõ¼þÊý Loadi=find(QD~=0);
%Loadi=[1:Busnum]';
RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1; %约束条件数
t_Bal_volt=Volt(Balance); t_Bal_volt=Volt(Balance);
Volt=sparse(ones(1,Busnum)); Volt=sparse(ones(1,Busnum));
Volt(Balance)=t_Bal_volt; Volt(Balance)=t_Bal_volt;
%Volt(Balance)=1;
UAngel=sparse(1,Busnum); UAngel=sparse(1,Busnum);
Init_Z=sparse(ones(1,RestraintCount)); Init_Z=sparse(ones(1,RestraintCount));
Init_W=sparse(-1*ones(1,RestraintCount)); Init_W=sparse(-1*ones(1,RestraintCount));
Init_L=sparse(ones(1,RestraintCount)); Init_L=sparse(ones(1,RestraintCount));
Init_U=sparse(ones(1,RestraintCount)); Init_U=sparse(ones(1,RestraintCount));
Init_Y=sparse(1,2*Busnum);%与学姐一致 Init_Y=sparse(1,2*Busnum);%与学姐一致
%Init_Y(1:2:2*Busnum)=1e-10;
%Init_Y(2:2:2*Busnum)=-1e-10;
tPU=sparse(GenU(:,2));% 发电机有功上限 tPU=sparse(GenU(:,2));% 发电机有功上限
tQU=sparse(PVQU(:,1));% 无功上限 tQU=sparse(PVQU(:,1));% 无功上限
tPL=sparse(GenL(:,2));% 发电机有功下限 tPL=sparse(GenL(:,2));% 发电机有功下限
tQL=sparse(PVQL(:,1));% 无功下限 tQL=sparse(PVQL(:,1));% 无功下限
%PG(4:5)=[4.5 4.5];
PG(PGi)=(tPU+tPL)/2; PG(PGi)=(tPU+tPL)/2;
%QG(4:5)=[0 1.45];
QG(PVi)=(tQU+tQL)/2; QG(PVi)=(tQU+tQL)/2;
wG=ones(size(PGi,1),1); wG=ones(size(PGi,1),1);
wD=ones(Busnum,1); randInt=randperm(size(Loadi,1));
%asc_randInt=sort(randInt,2,'ascend' );
%randPDind=randInt(1:200);
randPDind=randInt(1:0);
wD=ones(size(Loadi,1),1);
wD(randPDind)=0;%一些负荷不约束
%wD(Balance)=0; %wD(Balance)=0;
%PD0(randPDind)=0;
PD=1*PD0; PD=1*PD0;
%PD(PD==0)=.2; %PD(PD==0)=.2;

4
ObjectiveFun.m
View File

@ -1,6 +1,6 @@
function [out_arg]=ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD) function [out_arg]=ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
t1=PG(PGi)-PG0(PGi); t1=PG(PGi)-PG0(PGi);
t2=wG.*t1.*t1; t2=wG.*t1.*t1;
t3=wD.*((PD-PD0).*(PD-PD0)); t3=wD.*((PD(Loadi)-PD0(Loadi)).*(PD(Loadi)-PD0(Loadi)));
out_arg= sum(t2)+sum(t3); out_arg= sum(t2)+sum(t3);
end end

6
SolveIt.m
View File

@ -1,4 +1,4 @@
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) 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); 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;%+ddg*(Init_Z'+Init_W');
t1=diag(Init_L.\Init_Z-Init_U.\Init_W); t1=diag(Init_L.\Init_Z-Init_U.\Init_W);
@ -9,13 +9,13 @@ aa=[
]; ];
yy=[LxComa;-Ly]; yy=[LxComa;-Ly];
%% 平衡节点电压不变 %% 平衡节点电压不变
t=size(PVi,1)+size(PGi,1)+Busnum; t=size(PVi,1)+size(PGi,1)+size(Loadi,1);
aa(t+Balance,:)=0; aa(t+Balance,:)=0;
aa(:,t+Balance)=0; aa(:,t+Balance)=0;
aa(t+Balance,t+Balance)=1; aa(t+Balance,t+Balance)=1;
deltG(t+Balance,:)=0; deltG(t+Balance,:)=0;
%% %%
t=size(PVi,1)+size(PGi,1)+Busnum*2; t=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1;
aa(t+Balance,:)=0; aa(t+Balance,:)=0;
aa(:,t+Balance)=0; aa(:,t+Balance)=0;
aa(t+Balance,t+Balance)=1; aa(t+Balance,t+Balance)=1;

4
func_ddg.m
View File

@ -1,5 +1,5 @@
function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount) function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi)
ddg=sparse(size(PVi,1)+size(PGi,1)+2*Busnum,RestraintCount); ddg=sparse(size(PVi,1)+size(PGi,1)+size(Loadi,1)+1*Busnum,RestraintCount);
end end

9
func_ddh3.m
View File

@ -1,6 +1,6 @@
function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle) function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi)
%决定用循环重写 %决定用循环重写
ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
%AngleIJ=AngleIJMat-angle(GB); %AngleIJ=AngleIJMat-angle(GB);
mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
@ -47,8 +47,9 @@ t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV ;
]; ];
sizePGi=size(PGi,1); sizePGi=size(PGi,1);
sizePVi=size(PVi,1); sizePVi=size(PVi,1);
sizeLoadi=size(Loadi,1);
ddh=[ ddh=[
sparse(sizePGi+sizePVi+Busnum,ContrlCount); sparse(sizePGi+sizePVi+sizeLoadi,ContrlCount);
sparse(2*Busnum,sizePVi+sizePGi+Busnum),-t; sparse(2*Busnum,sizePVi+sizePGi+sizeLoadi),-t;
]; ];
end end

4
func_deltF.m
View File

@ -1,4 +1,4 @@
function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum) function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi)
%t1=PG(setdiff(PVi,Balance)); %t1=PG(setdiff(PVi,Balance));
% t2=Volt'*Volt; % t2=Volt'*Volt;
% t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); % t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
@ -8,7 +8,7 @@ function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum)
% PPG=([PQ(1),PBal])';%暂时用土办法处理一下 % PPG=([PQ(1),PBal])';%暂时用土办法处理一下
%% %%
t1=2*wG.*(PG(PGi) - PG0(PGi) ); t1=2*wG.*(PG(PGi) - PG0(PGi) );
t2=2*wD.*(PD-PD0); t2=2*wD.*(PD(Loadi)-PD0(Loadi));
deltF=[ deltF=[
sparse(t1); sparse(t1);
sparse(size(PVi,1),size(PVi,2)); sparse(size(PVi,1),size(PVi,2));

37
func_deltG.m
View File

@ -1,46 +1,27 @@
function deltG=func_deltG(Busnum,PVi,PGi) function deltG=func_deltG(Busnum,PVi,PGi,Loadi)
%% %%
sizePGi=size(PGi,1); sizePGi=size(PGi,1);
sizePVi=size(PVi,1); sizePVi=size(PVi,1);
sizeLoadi=size(Loadi,1);
%% %%
% dg1_dPg=eye(size(PGi,1));
% dg2_dPg=zeros(size(PGi,1),size(PVi,1));
% dg3_dPg=zeros(size(PGi,1),Busnum);
% dg4_dPg=zeros(size(PGi,1),Busnum);
dg1_dPg=sparse(1:sizePGi,1:sizePGi,ones(sizePGi,1),sizePGi,sizePGi); dg1_dPg=sparse(1:sizePGi,1:sizePGi,ones(sizePGi,1),sizePGi,sizePGi);
dg2_dPg=sparse(sizePGi,sizePVi); dg2_dPg=sparse(sizePGi,sizePVi);
dg3_dPg=sparse(sizePGi,Busnum); dg3_dPg=sparse(sizePGi,sizeLoadi);
dg4_dPg=sparse(sizePGi,Busnum); dg4_dPg=sparse(sizePGi,Busnum);
%% %%
% dg1_dQr=zeros(size(PVi,1),size(PGi,1));
% dg2_dQr=eye(size(PVi,1));
% dg3_dQr=zeros(size(PVi,1),Busnum);
% dg4_dQr=zeros(size(PVi,1),Busnum);
dg1_dQr=sparse(sizePVi,sizePGi); dg1_dQr=sparse(sizePVi,sizePGi);
dg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi); dg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi);
dg3_dQr=sparse(sizePVi,Busnum); dg3_dQr=sparse(sizePVi,sizeLoadi);
dg4_dQr=sparse(sizePVi,Busnum); dg4_dQr=sparse(sizePVi,Busnum);
%% %%
% dg1_dPD=zeros(Busnum,size(PGi,1)); dg1_dPD=sparse(size(Loadi,1),size(PGi,1));
% dg2_dPD=zeros(Busnum,size(PVi,1)); dg2_dPD=sparse(size(Loadi,1),size(PVi,1));
% dg3_dPD=eye(Busnum,Busnum); dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
% dg4_dPD=zeros(Busnum,Busnum); dg4_dPD=sparse(size(Loadi,1),Busnum);
dg1_dPD=sparse(Busnum,size(PGi,1));
dg2_dPD=sparse(Busnum,size(PVi,1));
dg3_dPD=sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
dg4_dPD=sparse(Busnum,Busnum);
%% %%
% dg1_dx=zeros(2*Busnum,size(PGi,1));
% dg2_dx=zeros(2*Busnum,size(PVi,1));
% dg3_dx=zeros(2*Busnum,Busnum);
% dg4_dx=zeros(2*Busnum,Busnum);
% for I=1:Busnum
% %dg3_dx(2*I,I)=1;ÔÝĘą¸ÄŇťĎÂ
% dg4_dx(I,I)=1;
% end
dg1_dx=sparse(2*Busnum,sizePGi); dg1_dx=sparse(2*Busnum,sizePGi);
dg2_dx=sparse(2*Busnum,sizePVi); dg2_dx=sparse(2*Busnum,sizePVi);
dg3_dx=sparse(2*Busnum,Busnum); dg3_dx=sparse(2*Busnum,sizeLoadi);
dg4_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum); dg4_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
sparse(Busnum,Busnum); sparse(Busnum,Busnum);
]; ];

22
func_deltH.m
View File

@ -1,24 +1,8 @@
function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle) function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle,Loadi)
% dH_dPg=zeros(size(PGi,1),2*Busnum);
%
% 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_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum); dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum);
% dH_dQr=zeros(size(PVi,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)+Busnum)=1;
% end
dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,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:Busnum,1:Busnum,-ones(Busnum,1),Busnum,Busnum) sparse(Busnum,Busnum)]; dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)];
%Angle=angle(GB); dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c); %形成雅克比矩阵
dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat,UAngel,r,c); %形成雅克比矩阵
%deltH=[dH_dPg;dH_dQr;dH_dx'];%dH_dx 需要使用一下转置 暂时改一下
deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx']; deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx'];
end end

2
func_deltdeltF.m
View File

@ -1,6 +1,6 @@
function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0) function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0)
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']; C=[wG' zeros(size(PVi))' wD'];
sizeC=size(C,2); sizeC=size(C,2);
diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC); diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC);

2
jacobian_M3.m
View File

@ -1,4 +1,4 @@
function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat,UAngel,r,c) function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c)
%************************************************************************** %**************************************************************************
% 程序功能 : 子函数——形成雅可比矩阵Jacobian % 程序功能 : 子函数——形成雅可比矩阵Jacobian
% 编 者： % 编 者：