diff --git a/Run_YALMIP.asv b/Run_YALMIP.asv index 628fdb5..1f8e64f 100644 --- a/Run_YALMIP.asv +++ b/Run_YALMIP.asv @@ -5,7 +5,7 @@ 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, ... Liner,Linex,Lineb,Transforr,Transforx,Transfork0]= ... - pf('ieee30.dat'); + pf('E:/算例/柳金Ⅰ926_21671693_2012-09-06/newFIle20.txt'); %% 潮流等式 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); @@ -24,48 +24,51 @@ QG0(Balance)=QGBal(Balance); QG0(PVi)=QGBal(PVi); PG(Balance)=PGBal(Balance); QG(PVi)=QGBal(PVi); -PF=1; -AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); -dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt'; -dPD=abs(dP./PD); -dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt'; -dQD=abs(dQ./QD); -maxdPQ=max([dPD(dPD<10);dQD(dQD<10)]); - -[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD); - +Volt0=Volt; +% PF=1; +% AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); +% dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt'; +% dPD=abs(dP./PD); +% dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt'; +% dQD=abs(dQ./QD); +% maxdPQ=max([dPD(dPD<10);dQD(dQD<10)]); +[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD); %% 定义变量 -Volt=sdpvar(Busnum,1); -UAngel=sdpvar(Busnum,1); +BalVolt=Volt(Balance); +% Volt=sdpvar(Busnum,1); +% UAngel=sdpvar(Busnum,1); % PG=sdpvar(Busnum,1); % QG=sdpvar(Busnum,1); PD=sdpvar(Busnum,1); QD=sdpvar(Busnum,1); AngleIJ=sdpvar(Busnum,Busnum,'full'); +%% 加噪音 +PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0; +QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0; %% 目标函数 -Objective=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi); +Objective=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,Volt,Volt0,wPG,wQG,wPD,wQD,wVolt,Loadi); %AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); %% 赋初值,可以加快求解速度。 assign(Volt(:),1); assign(UAngel(:),0); assign(PD(:),PD0(:)); assign(QD(:),QD0(:)); -% assign(PG(:),PG0(:)); -% assign(QG(:),QG0(:)); %% YALMIP部分 dP=PG0-PD-diag(Volt)*Y.*cos( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; dQ=QG0-QD-diag(Volt)*Y.*sin( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; -Loadi=PD0~=0 | QD0~=0; +Loadi=PD0~=0 | QD0~=0 |PG0~=0|QG0~=0; Constraints = [%AngleIJ-sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum)==0, ... - dP==0, ... - dQ==0, ... - PG(setxor(1:Busnum, PGi) )==0, ... - QG(setxor(1:Busnum, PVi) )==0, ... + dP(setdiff(1:Busnum,Loadi))==0, ... + dQ(setdiff(1:Busnum,Loadi))==0, ... +% dP==0, ... +% dQ==0, ... PD(PD0==0)==0, ... QD(QD0==0)==0, ... 0.9*ones(Busnum,1)<=Volt<=1.1*ones(Busnum,1), ... - -60*ones(Busnum,1)<=PG<=60*ones(Busnum,1), ... - -60*ones(Busnum,1)<=QG<=60*ones(Busnum,1) + Volt(Balance)==BalVolt, ... + UAngel(Balance)==0, ... + 0.8*PD0<=PD<=1.2*PD0; + 0.8*QD0<=QD<=1.2*QD0; ]; options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt','usex0','1'); sol = solvesdp(Constraints,Objective,options); @@ -76,7 +79,6 @@ if sol.problem == 0 dVangle=double(UAngel) fprintf('ojb\n'); optimalObj=double(Objective) - double(PG)-PG0; sol else display('Hmm, something went wrong!'); diff --git a/Run_YALMIP.m b/Run_YALMIP.m index b098f5b..710cd91 100644 --- a/Run_YALMIP.m +++ b/Run_YALMIP.m @@ -35,8 +35,8 @@ Volt0=Volt; [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD); %% 定义变量 BalVolt=Volt(Balance); -Volt=sdpvar(Busnum,1); -UAngel=sdpvar(Busnum,1); +% Volt=sdpvar(Busnum,1); +% UAngel=sdpvar(Busnum,1); % PG=sdpvar(Busnum,1); % QG=sdpvar(Busnum,1); PD=sdpvar(Busnum,1); @@ -46,45 +46,45 @@ AngleIJ=sdpvar(Busnum,Busnum,'full'); PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0; QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0; %% 目标函数 -Objective=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,Volt,Volt0,wPG,wQG,wPD,wQD,wVolt,Loadi); +% Objective=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,Volt,Volt0,wPG,wQG,wPD,wQD,wVolt,Loadi); %AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); %% 赋初值,可以加快求解速度。 -assign(Volt(:),1); -assign(UAngel(:),0); -assign(PD(:),PD0(:)); -assign(QD(:),QD0(:)); -%% YALMIP部分 -dP=PG0-PD-diag(Volt)*Y.*cos( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; -dQ=QG0-QD-diag(Volt)*Y.*sin( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; -Loadi=PD0~=0 | QD0~=0 |PG0~=0|QG0~=0; -Constraints = [%AngleIJ-sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum)==0, ... - dP(setdiff(1:Busnum,Loadi))==0, ... - dQ(setdiff(1:Busnum,Loadi))==0, ... -% dP==0, ... -% dQ==0, ... - PD(PD0==0)==0, ... - QD(QD0==0)==0, ... - 0.9*ones(Busnum,1)<=Volt<=1.1*ones(Busnum,1), ... - Volt(Balance)==BalVolt, ... - UAngel(Balance)==0, ... - 0.8*PD0<=PD<=1.2*PD0; - 0.8*QD0<=QD<=1.2*QD0; - ]; -options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt','usex0','1'); -sol = solvesdp(Constraints,Objective,options); -if sol.problem == 0 - fprintf('Volt\n'); - dvolt=double(Volt) - fprintf('VoltAngle\n'); - dVangle=double(UAngel) - fprintf('ojb\n'); - optimalObj=double(Objective) - sol -else - display('Hmm, something went wrong!'); - sol.info - sol.solveroutput - yalmiperror(sol.problem) -end +% assign(Volt(:),1); +% assign(UAngel(:),0); +% assign(PD(:),PD0(:)); +% assign(QD(:),QD0(:)); +% %% YALMIP部分 +% dP=PG0-PD-diag(Volt)*Y.*cos( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; +% dQ=QG0-QD-diag(Volt)*Y.*sin( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; +% Loadi=PD0~=0 | QD0~=0 |PG0~=0|QG0~=0; +% Constraints = [%AngleIJ-sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum)==0, ... +% dP(setdiff(1:Busnum,Loadi))==0, ... +% dQ(setdiff(1:Busnum,Loadi))==0, ... +% % dP==0, ... +% % dQ==0, ... +% PD(PD0==0)==0, ... +% QD(QD0==0)==0, ... +% 0.9*ones(Busnum,1)<=Volt<=1.1*ones(Busnum,1), ... +% Volt(Balance)==BalVolt, ... +% UAngel(Balance)==0, ... +% 0.8*PD0<=PD<=1.2*PD0; +% 0.8*QD0<=QD<=1.2*QD0; +% ]; +% options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt','usex0','1'); +% sol = solvesdp(Constraints,Objective,options); +% if sol.problem == 0 +% fprintf('Volt\n'); +% dvolt=double(Volt) +% fprintf('VoltAngle\n'); +% dVangle=double(UAngel) +% fprintf('ojb\n'); +% optimalObj=double(Objective) +% sol +% else +% display('Hmm, something went wrong!'); +% sol.info +% sol.solveroutput +% yalmiperror(sol.problem) +% end toc