From 45f7dc7dff05037442e17ec0e44dbdddd8ef7cae Mon Sep 17 00:00:00 2001 From: facat Date: Wed, 17 Apr 2013 10:04:43 +0800 Subject: [PATCH] =?UTF-8?q?=E4=BB=A3=E7=A0=81=E5=86=99=E5=A5=BD=E4=BA=86?= =?UTF-8?q?=EF=BC=8C=E6=AD=A3=E5=9C=A8=E8=B0=83=E8=AF=95?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- @Opti/Geteb.m | 7 ++++ @Opti/Getgle.m | 7 ++++ @Opti/Opti.m | 31 ++++++++++++++ @Opti/equ.asv | 27 +++++++++++++ @Opti/equ.m | 30 ++++++++++++++ @Opti/init.m | 19 +++++++++ @Opti/obj.asv | 15 +++++++ @Opti/obj.m | 17 ++++++++ OPF_Init.m | 6 +-- Run_YALMIP.asv | 107 ++++++++++++++++++++++++++++++------------------- Run_YALMIP.m | 32 +++++++++++++-- 11 files changed, 251 insertions(+), 47 deletions(-) create mode 100644 @Opti/Geteb.m create mode 100644 @Opti/Getgle.m create mode 100644 @Opti/Opti.m create mode 100644 @Opti/equ.asv create mode 100644 @Opti/equ.m create mode 100644 @Opti/init.m create mode 100644 @Opti/obj.asv create mode 100644 @Opti/obj.m diff --git a/@Opti/Geteb.m b/@Opti/Geteb.m new file mode 100644 index 0000000..244762a --- /dev/null +++ b/@Opti/Geteb.m @@ -0,0 +1,7 @@ +function [ output_args ] = Geteb( this ) +%BE Summary of this function goes here +% Detailed explanation goes here +output_args=this.eb; + +end + diff --git a/@Opti/Getgle.m b/@Opti/Getgle.m new file mode 100644 index 0000000..f12771c --- /dev/null +++ b/@Opti/Getgle.m @@ -0,0 +1,7 @@ +function [ output_args ] = Getgle( this ) +%GLE Summary of this function goes here +% Detailed explanation goes here +output_args=this.gle; + +end + diff --git a/@Opti/Opti.m b/@Opti/Opti.m new file mode 100644 index 0000000..ed0fa1a --- /dev/null +++ b/@Opti/Opti.m @@ -0,0 +1,31 @@ +classdef Opti + %OPTI Summary of this class goes here + % Detailed explanation goes here + + properties + PD0=NaN; + QD0=NaN; + Volt0=NaN; + wPD=NaN; + wQD=NaN; + wVolt=NaN; + PDi=NaN; + QDi=NaN; + Y=NaN; + Angel=NaN; + r=NaN; + c=NaN; + PG=NaN; + QG=NaN; + gle=NaN;%opti toolbox大于小于等于 + eb=NaN;%等式、不等式约束右边的值 + Balance=NaN; + end + + methods + function this=Opti() + end + end + +end + diff --git a/@Opti/equ.asv b/@Opti/equ.asv new file mode 100644 index 0000000..7abd80b --- /dev/null +++ b/@Opti/equ.asv @@ -0,0 +1,27 @@ +function [ output_args,this ] = equ( this ) +%EQU Summary of this function goes here +% Detailed explanation goes here +PDi=this.PDi; +QDi=this.QDi; +Volt0=this.Volt0; +Angel=this.Angel; +r=this.r; +c=this.c; +PG=this.PG; +QG=this.QG; +PD=x(1:length(PDi)); +QD=x(length(PDi)+1:length(PDi)+length(QDi)); +Volt=x(length(PDi)+length(QDi)+1:length(PDi)+length(QDi)+length(Volt0)); +VAngel=x(length(PDi)+length(QDi)+length(Volt0):end); +PD_=zeros(PD,length(Volt0)); +QD_=zeros(QD,length(Volt0)); +PD_(PDi)=PD; +QD_(QDi)=QD; +busNum=length(Volt0); +VMatrix=sparse(r,c,VAngel(r)-VAngel(c)-Angel,busNum,busNum); +dP=PG-PD_-diag(Volt)*(Y.*cos(VMatrix))*Volt; +dQ=QG-QD_-diag(Volt)*(Y.*spfun(@sin,VMatrix))*Volt; +output_args=[dP;dQ]; +this.gle=ze +end + diff --git a/@Opti/equ.m b/@Opti/equ.m new file mode 100644 index 0000000..f994058 --- /dev/null +++ b/@Opti/equ.m @@ -0,0 +1,30 @@ +function [ output_args,this ] = equ( this,x ) +%EQU Summary of this function goes here +% Detailed explanation goes here +PDi=this.PDi; +QDi=this.QDi; +Volt0=this.Volt0; +Angel=this.Angel; +r=this.r; +c=this.c; +PG=this.PG; +QG=this.QG; +Balance=this.Balance; +Y=this.Y; +PD=x(1:length(PDi)); +QD=x(length(PDi)+1:length(PDi)+length(QDi)); +Volt=x(length(PDi)+length(QDi)+1:length(PDi)+length(QDi)+length(Volt0)); +VAngel=x(length(PDi)+length(QDi)+length(Volt0):end); +PD_=zeros(length(Volt0),1); +QD_=zeros(length(Volt0),1); +PD_(PDi)=PD; +QD_(QDi)=QD; +busNum=length(Volt0); +VMatrix=sparse(r,c,VAngel(r)-VAngel(c)-Angel,busNum,busNum); +dP=PG-PD_-diag(Volt)*(Y.*cos(VMatrix))*Volt; +dQ=QG-QD_-diag(Volt)*(Y.*spfun(@sin,VMatrix))*Volt; +output_args=[dP;dQ]; +this.gle=zeros(length(output_args),1); +this.eb=this.gle; +end + diff --git a/@Opti/init.m b/@Opti/init.m new file mode 100644 index 0000000..a36e034 --- /dev/null +++ b/@Opti/init.m @@ -0,0 +1,19 @@ +function [ this ] = init( this,Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angel,r,c,PG,QG,Balance ) +%INIT Summary of this function goes here +% Detailed explanation goes here +this.PDi=PDi; +this.QDi=QDi; +this.wPD=wPD; +this.wQD=wQD; +this.PD0=PD0; +this.QD0=QD0; +this.Y=Y; +this.Angel=Angel; +this.r=r; +this.c=c; +this.PG=PG; +this.QG=QG; +this.Balance=Balance; +this.Volt0=Volt0; +end + diff --git a/@Opti/obj.asv b/@Opti/obj.asv new file mode 100644 index 0000000..a1efcc8 --- /dev/null +++ b/@Opti/obj.asv @@ -0,0 +1,15 @@ +function [ output_args ] = obj( this ) +%OBJ Summary of this function goes here +% Detailed explanation goes here +wPD=this.wPD; +wQD=this.wQD; +PD=this.PD; +QD=this.QD; +wVolt=t +t4=wPD.*(PD(PD0~=0)-PD0(PD0~=0)).^2; +t5=wQD.*(QD(QD0~=0)-QD0(QD0~=0)).^2; +t6=wVolt.*(Volt-Volt0').^2; +out_arg=sum(t4)+sum(t5)+sum(t6); + +end + diff --git a/@Opti/obj.m b/@Opti/obj.m new file mode 100644 index 0000000..4811ef5 --- /dev/null +++ b/@Opti/obj.m @@ -0,0 +1,17 @@ +function [ output_args ] = obj( this,x ) +%OBJ Summary of this function goes here +% Detailed explanation goes here +PDi=this.PDi; +QDi=this.QDi; +PD=x(1:length(PDi)); +QD=x(length(PDi)+1:length(PDi)+length(QDi)); +wPD=this.wPD; +wQD=this.wQD; +% wVolt=this.wVolt; +PD0=this.PD0; +QD0=this.QD0; +t4=(wPD(PDi).*(PD-PD0(PDi))).^2; +t5=(wQD(QDi).*(QD-QD0(QDi))).^2; +output_args=sum(t4)+sum(t5); +end + diff --git a/OPF_Init.m b/OPF_Init.m index 3c4db83..91efd95 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -33,9 +33,9 @@ wQG=1*ones(size(PVi,1),1); %randInt=randperm(size(Loadi,1)); %randPDind=randInt(1:10); randPDind=0; -wPD=1/0.05^2*ones(sum(PD~=0),1); -wQD=1/0.05^2*ones(sum(QD~=0),1); -wVolt=1/0.05^2*ones(Busnum,1); +wPD=1/0.05*ones(length(PD),1); +wQD=1/0.05*ones(length(QD),1); +wVolt=1/0.05*ones(Busnum,1); %wD(randPDind)=0;%一些负荷不约束 %wD(7)=0; % wD(11)=0; diff --git a/Run_YALMIP.asv b/Run_YALMIP.asv index 1f8e64f..0096424 100644 --- a/Run_YALMIP.asv +++ b/Run_YALMIP.asv @@ -39,52 +39,77 @@ BalVolt=Volt(Balance); % UAngel=sdpvar(Busnum,1); % PG=sdpvar(Busnum,1); % QG=sdpvar(Busnum,1); -PD=sdpvar(Busnum,1); -QD=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,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 +%% Opti Toolbox +Busnum=length(Volt); +PDi=find(PD~=0); +QDi=find(QD~=0); +% PD0=PD0(PDi); +% QD0=QD0(QDi); +seOpti=Opti(); +seOpti=seOpti.init(PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle); +opts = optiset('solver','filtersd'); +opts.maxiter=85500; +opts.maxtime=3000; +opts.maxfeval=85000; +opts.maxnodes=85000; +opts.tolrfun=1e-4; +opts.tolafun=1e-4; +opts.warnings='all'; +opts.display='off'; +x0=zeros(length(PDi)+length(QDi),1); +% [~,seOpti]=seOpti.equ(x0); +% nlrhs=seOpti.nlrhs(); +% nle=seOpti.nle(); +% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) +Opt = opti('fun',@seOpti.obj,'ndec',length(PDi)+length(QDi),'options',opts) +[x,fval,exitflag,info] = solve(Opt,x0); +info +fprintf('目标函数: %f\n',fval); toc diff --git a/Run_YALMIP.m b/Run_YALMIP.m index 710cd91..193d57e 100644 --- a/Run_YALMIP.m +++ b/Run_YALMIP.m @@ -39,8 +39,8 @@ BalVolt=Volt(Balance); % UAngel=sdpvar(Busnum,1); % PG=sdpvar(Busnum,1); % QG=sdpvar(Busnum,1); -PD=sdpvar(Busnum,1); -QD=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; @@ -86,5 +86,31 @@ QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0; % sol.solveroutput % yalmiperror(sol.problem) % end - +%% Opti Toolbox +Busnum=length(Volt); +PDi=find(PD~=0); +QDi=find(QD~=0); +% PD0=PD0(PDi); +% QD0=QD0(QDi); +seOpti=Opti(); +seOpti=seOpti.init(Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle,r,c,PG,QG,Balance); +opts = optiset('solver','filtersd'); +opts.maxiter=85500; +opts.maxtime=3000; +opts.maxfeval=85000; +opts.maxnodes=85000; +opts.tolrfun=1e-4; +opts.tolafun=1e-4; +opts.warnings='all'; +opts.display='off'; +x0=zeros(length(PDi)+length(QDi)+length(Volt)*2,1); +% [~,seOpti]=seOpti.equ(x0); +[~,seOpti]=seOpti.equ(x0); +nlrhs=seOpti.Geteb(); +nle=seOpti.Getgle(); +Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2+length(PDi)+length(QDi),'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) +% Opt = opti('fun',@seOpti.obj,'ndec',length(PDi)+length(QDi),'options',opts) +[x,fval,exitflag,info] = solve(Opt,x0); +info +fprintf('目标函数: %f\n',fval); toc