From b7498c2b15967efa8088b18c03e1cccad915e506 Mon Sep 17 00:00:00 2001 From: "dmy@lab" Date: Thu, 19 Mar 2015 13:19:57 +0800 Subject: [PATCH] =?UTF-8?q?=E6=94=B9=E6=88=90=E4=BA=86=E5=90=AB0=E6=B3=A8?= =?UTF-8?q?=E5=85=A5=E7=BA=A6=E6=9D=9F=EF=BC=88=E4=B8=8D=E6=98=AF=E6=BD=AE?= =?UTF-8?q?=E6=B5=81=E7=BA=A6=E6=9D=9F=EF=BC=89=E7=9A=84=E6=9C=80=E5=B0=8F?= =?UTF-8?q?=E4=BA=8C=E4=B9=98=E6=B3=95=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dmy@lab --- @Opti/equ.asv | 41 ++------------------------- @Opti/equ.m | 78 +++++++++++++-------------------------------------- @Opti/obj.asv | 41 ++++++++++++++++----------- @Opti/obj.m | 35 +++++++++++++++-------- OPF_Init.m | 7 +++-- Run_YALMIP2.m | 16 +++++------ 6 files changed, 83 insertions(+), 135 deletions(-) diff --git a/@Opti/equ.asv b/@Opti/equ.asv index a09d469..d1f2b87 100644 --- a/@Opti/equ.asv +++ b/@Opti/equ.asv @@ -14,16 +14,13 @@ 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)+1:end); +Volt=x(1:length(Volt0)); +VAngel=x(length(Volt0)+1: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];%潮流约束 %% 使用零注入begin % zeroInj=setdiff(1:busNum,union(Balance,union(PDi,QDi))); @@ -35,38 +32,6 @@ output_args=[output_args;Volt(Balance)-1]; output_args=[output_args;VAngel(Balance)]; this.cu=zeros(length(output_args),1); this.cl=zeros(length(output_args),1); -%% 开始增加不等式约束-电压\PD\QD -%电压不等式约束 -output_args=[output_args;Volt]; -% this.cu=[this.cu;1.07*ones(length(Volt),1)];%9节点 -% this.cl=[this.cl;0.93*ones(length(Volt),1)];%9节点 -this.cu=[this.cu;1.00*ones(length(Volt),1)]; -this.cl=[this.cl;0.90*ones(length(Volt),1)]; -%% PD -rPD=this.rPD; -output_args=[output_args;PD]; -% PDU=[0.124;0.315;0.5;1;1;0.5;0.63;0.4];%原始数据 -% PDU=[0.63;0.4;0.5;1;0.8;1;0.63;0.4];%偏差比较大 -%Generate values from the uniform distribution on the interval [a, b]. -global loadFlag; -if loadFlag==0 - r =-0.201 + (-0.01-(-0.2)).*rand(length(rPD),1); -else - -end - -this.cu=[this.cu;(r+1.4).*rPD]; -this.cl=[this.cl;(r+1).*rPD]; -% this.cu=[this.cu;PDU]; -% this.cl=[this.cl;0*PDU]; -%% QD -rQD=this.rQD; -output_args=[output_args;QD]; -% QDU=PDU; -this.cu=[this.cu;(r+1.4).*rQD]; -this.cl=[this.cl;(r+1).*rQD]; -% this.cu=[this.cu;QDU]; -% this.cl=[this.cl;0*QDU]; %% 稠密化 output_args=full(output_args); this.cu=full(this.cu); diff --git a/@Opti/equ.m b/@Opti/equ.m index 2087110..83858f1 100644 --- a/@Opti/equ.m +++ b/@Opti/equ.m @@ -2,74 +2,36 @@ function [ output_args,this ] = equ( this,x ) %EQU Summary of this function goes here % Detailed explanation goes here output_args=[]; -PDi=this.PDi; -QDi=this.QDi; Volt0=this.Volt0; +Balance=this.Balance; +% Volt=x(1:length(Volt0)); +% VAngel=x(length(Volt0)+1:end); +%% 使用零注入begin +SEVolt=x(1:length(Volt0)); +SEVAngel=x(length(SEVolt)+1:end); 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)+1: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];%潮流约束 -%% 使用零注入begin -% zeroInj=setdiff(1:busNum,union(Balance,union(PDi,QDi))); -% output_args=[output_args(zeroInj);output_args(length(Volt0)+zeroInj)]; -%% 使用零注入end -%电压恒定 -output_args=[output_args;Volt(Balance)-1]; -%相角恒定 -output_args=[output_args;VAngel(Balance)]; +VMatrix=sparse(r,c,SEVAngel(r)-SEVAngel(c)-Angel,busNum,busNum); +zerosP=diag(SEVolt)*(Y.*cos(VMatrix))*SEVolt; +zerosQ=diag(SEVolt)*(Y.*spfun(@sin,VMatrix))*SEVolt; +PDi=this.PDi; +QDi=this.QDi; +zeroInj=setdiff(1:busNum,union(Balance,union(PDi,QDi))); +output_args=[zerosP(zeroInj);zerosQ(zeroInj)]; +this.cu=zeros(length(output_args),1); +this.cl=zeros(length(output_args),1); +%% 使用零注入end +%电压恒定 +output_args=[output_args;SEVolt(Balance)-1]; +%相角恒定 +output_args=[output_args;SEVAngel(Balance)]; this.cu=zeros(length(output_args),1); this.cl=zeros(length(output_args),1); -%% 开始增加不等式约束-电压\PD\QD -%电压不等式约束 -output_args=[output_args;Volt]; -% this.cu=[this.cu;1.07*ones(length(Volt),1)];%9节点 -% this.cl=[this.cl;0.93*ones(length(Volt),1)];%9节点 -this.cu=[this.cu;1.00*ones(length(Volt),1)]; -this.cl=[this.cl;0.90*ones(length(Volt),1)]; -%% PD -rPD=this.rPD; -output_args=[output_args;PD]; -% PDU=[0.124;0.315;0.5;1;1;0.5;0.63;0.4];%原始数据 -% PDU=[0.63;0.4;0.5;1;0.8;1;0.63;0.4];%偏差比较大 -%Generate values from the uniform distribution on the interval [a, b]. -global loadFlag; -if loadFlag==1 - r=load('rLD'); - r=r.r; - -else - %r=-0.2; - r =-0.201 + (-0.001-(-0.2)).*rand(length(rPD),1); -end - -this.cu=[this.cu;(r+1.4).*rPD]; -this.cl=[this.cl;(r+1).*rPD]; -% this.cu=[this.cu;PDU]; -% this.cl=[this.cl;0*PDU]; -%% QD -rQD=this.rQD; -output_args=[output_args;QD]; -% QDU=PDU; -this.cu=[this.cu;(r+1.4).*rQD]; -this.cl=[this.cl;(r+1).*rQD]; -% this.cu=[this.cu;QDU]; -% this.cl=[this.cl;0*QDU]; %% 稠密化 output_args=full(output_args); this.cu=full(this.cu); diff --git a/@Opti/obj.asv b/@Opti/obj.asv index 24ca46b..2ba4ee1 100644 --- a/@Opti/obj.asv +++ b/@Opti/obj.asv @@ -3,32 +3,39 @@ function [ output_args ] = obj( this,x ) % 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; PD0=this.PD0; QD0=this.QD0; -t4=(wPD(PDi).*(PD-PD0(PDi))).^2; -% t4=wPD(PDi).*((PD-PD0(PDi))./PD0(PDi)).^2; -t5=(wQD(QDi).*(QD-QD0(QDi))).^2; -% t5=wQD(QDi).*((QD-QD0(QDi))./QD0(QDi)).^2; Volt0=this.Volt0; -SEVolt=x(length(PDi)+length(QDi)+1:length(PDi)+length(QDi)+length(Volt0)); +SEVolt=x(1:length(Volt0)); wVolt=this.wVolt; -t6=(wVolt.*(SEVolt-Volt0)).^2; -% t6=wVolt.*((SEVolt-Volt0)./Volt0).^2; +t6=(wVolt.*((SEVolt-Volt0))).^2; Balance=this.Balance; t6(Balance)=0;%不计算平衡节点的 -t6(2:0) % 电流 -SEVAngel=x(length(PDi)+length(QDi)+length(SEVolt)+1:end); -loadCurrent=LoadCurrent( SEVolt,SEVAngel,PD0(PDi),QD0(QDi),PDi,QDi ); -wLoadCurrent=this.wLoadCurrent; -mLoadCurrent=this.mLoadCurrent; -t7=(wLoadCurrent.*(loadCurrent-mLoadCurrent)).^2; -t7=wLoadCurrent.*((loadCurrent-mLoadCurrent)./mLoadCurrent).^2; -output_args=sum(t4)+sum(t5)+sum(t6)+sum(0); +SEVAngel=x(length(SEVolt)+1:end); +busNum=length(Volt0); +Angel=this.Angel; +r=this.r; +c=this.c; +PG=this.PG; +QG=this.QG; + +VMatrix=sparse(r,c,SEVAngel(r)-SEVAngel(c)-Angel,busNum,busNum); +PD_=PG-diag(SEVolt)*(Y.*cos(VMatrix))*Volt; +QD_=QG-diag(SEVolt)*(Y.*spfun(@sin,VMatrix))*Volt; +t4=(wPD(PDi).*(PD_-PD0(PDi))).^2; +t5=(wQD(QDi).*(QD_-QD0(QDi))).^2; +% loadCurrent=LoadCurrent( SEVolt,SEVAngel,PD0(PDi),QD0(QDi),PDi,QDi ); +% wLoadCurrent=this.wLoadCurrent; +% mLoadCurrent=this.mLoadCurrent; +% t7=(wLoadCurrent.*(loadCurrent-mLoadCurrent)).^2; +% t7=wLoadCurrent.*((loadCurrent-mLoadCurrent)./mLoadCurrent).^2; +output_args=sum(t4)+sum(t5)+sum(t6); output_args=full(output_args); + + + end diff --git a/@Opti/obj.m b/@Opti/obj.m index 3e971e9..1698e78 100644 --- a/@Opti/obj.m +++ b/@Opti/obj.m @@ -3,28 +3,41 @@ function [ output_args ] = obj( this,x ) % 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; PD0=this.PD0; QD0=this.QD0; -t4=(wPD(PDi).*(PD-PD0(PDi))).^2; -t5=(wQD(QDi).*(QD-QD0(QDi))).^2; Volt0=this.Volt0; -SEVolt=x(length(PDi)+length(QDi)+1:length(PDi)+length(QDi)+length(Volt0)); +SEVolt=x(1:length(Volt0)); wVolt=this.wVolt; +% wVolt=ones(length(wVolt),1); t6=(wVolt.*((SEVolt-Volt0))).^2; Balance=this.Balance; t6(Balance)=0;%不计算平衡节点的 % 电流 -SEVAngel=x(length(PDi)+length(QDi)+length(SEVolt)+1:end); -loadCurrent=LoadCurrent( SEVolt,SEVAngel,PD0(PDi),QD0(QDi),PDi,QDi ); -wLoadCurrent=this.wLoadCurrent; -mLoadCurrent=this.mLoadCurrent; +SEVAngel=x(length(SEVolt)+1:end); +busNum=length(Volt0); +Angel=this.Angel; +r=this.r; +c=this.c; +PG=this.PG; +QG=this.QG; +Y=this.Y; +VMatrix=sparse(r,c,SEVAngel(r)-SEVAngel(c)-Angel,busNum,busNum); +PD_=PG-diag(SEVolt)*(Y.*cos(VMatrix))*SEVolt; +QD_=QG-diag(SEVolt)*(Y.*spfun(@sin,VMatrix))*SEVolt; +t4=(wPD(PDi).*(PD_(PDi)-PD0(PDi))).^2; +t5=(wQD(QDi).*(QD_(QDi)-QD0(QDi))).^2; +% loadCurrent=LoadCurrent( SEVolt,SEVAngel,PD0(PDi),QD0(QDi),PDi,QDi ); +% wLoadCurrent=this.wLoadCurrent; +% mLoadCurrent=this.mLoadCurrent; % t7=(wLoadCurrent.*(loadCurrent-mLoadCurrent)).^2; -t7=wLoadCurrent.*((loadCurrent-mLoadCurrent)./mLoadCurrent).^2; -output_args=sum(t4)+sum(t5)+sum(t6); +% t7=wLoadCurrent.*((loadCurrent-mLoadCurrent)./mLoadCurrent).^2; +% output_args=sum(t4)+sum(t5)+sum(t6); +output_args=sum(t6)+sum(t4)+sum(t5); output_args=full(output_args); + + + end diff --git a/OPF_Init.m b/OPF_Init.m index bf57e09..9a97f3d 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -5,6 +5,7 @@ PDi=find(rPD~=0); %QDi=find(rQD~=0); QDi=PDi; notLoadi=setdiff(1:Busnum,Loadi); +notLoadi=[]; PGOnly=setdiff(PGi,PDi); QGOnly=setdiff(PVi,QDi); PDOnly=setdiff(PDi,PGi); @@ -107,9 +108,9 @@ wVolt(notLoadi)=0; % wVolt(noLoad)=0; %% 0% noLoad=[]; -wPD(1:end)=0; -wQD(1:end)=0; -wVolt(1:end)=0; +% wPD(1:end)=0; +% wQD(1:end)=0; +% wVolt(1:end)=0; %wD(randPDind)=0;%一些负荷不约束 %wD(7)=0; % wD(11)=0; diff --git a/Run_YALMIP2.m b/Run_YALMIP2.m index 77bcaa1..7486365 100644 --- a/Run_YALMIP2.m +++ b/Run_YALMIP2.m @@ -9,7 +9,8 @@ for I=1:1 close all; [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('E:\算例\feeder33\feeder33.txt'); + Liner,Linex,Lineb,Transforr,Transforx,Transfork0]=pf('E:/算例/东际911_2751267_2012-09-05/newFIle20-使用.txt'); +% 'E:\算例\feeder33\feeder33.txt' %pf('C:\bpa\长虹世纪线_556844_2013-10-28\长虹世纪线_556844_2013-10-28_iPso_newFil %e.txt'); % pf('E:/算例/东际911_2751267_2012-09-05/pyth-增加3倍.txt'); 20131125 @@ -93,14 +94,13 @@ for I=1:1 opts.tolafun=1e-4; opts.warnings='all'; opts.display='off'; - x0=[rand()*rPD(PDi);rand()*rQD(QDi); ... - ones(length(Volt),1); ... + x0=[ones(length(Volt),1); ... zeros(length(Volt),1)]; % x0=[PD(PDi);QD(QDi);xVolt';xUAngel']; [~,seOpti]=seOpti.equ(x0); cl=seOpti.Getcl(); cu=seOpti.Getcu(); - Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2+length(PDi)+length(QDi),'nl',@seOpti.equ,cl,cu,'options',opts) + Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2,'nl',@seOpti.equ,cl,cu,'options',opts) % Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2+length(PDi)+length(QDi),'options',opts) [x,fval,exitflag,info] = solve(Opt,x0); info @@ -109,8 +109,8 @@ for I=1:1 toc rVolt=Volt0'; rVAngel=xUAngel'; - SEVolt=x(length(PDi)+length(QDi)+1:length(Volt)+length(PDi)+length(QDi)); - SEVAngel=x(length(Volt)+length(PDi)+length(QDi)+1:end); + SEVolt=x(1:length(Volt)); + SEVAngel=x(length(Volt)+1:end); fprintf('最大偏差\n') PD=x(1:length(PDi)); QD=x(length(PDi)+1:length(PDi)+length(QDi)); @@ -133,9 +133,9 @@ for I=1:1 fprintf('真实值统计误差为: %f',stE); stE=StErrorZ(SEVolt,PD,QD,mVolt,mPD(PDi),mQD(QDi),noLoadi); fprintf('测量值统计误差为: %f',stE); - NormalizedResiduals(x,sigma,PDi,QDi,Volt0,mPD,mQD); +% NormalizedResiduals(x,sigma,PDi,QDi,Volt0,mPD,mQD); SECurrent=LoadCurrent( SEVolt,SEVAngel,PD,QD,PDi,QDi ); - [flag,t1,t2]=MaxSigma( x,PDi,QDi,mPD,mQD,mVolt,mLoadCurrent,sigma,rPD,rQD,rVolt,rLoadCurrent ); +% [flag,t1,t2]=MaxSigma( x,PDi,QDi,mPD,mQD,mVolt,mLoadCurrent,sigma,rPD,rQD,rVolt,rLoadCurrent ); figure(); SEMeasDeivation(PD,full(PD0(PDi)),QD,full(QD0(QDi)),rPD(PDi),rQD(QDi)); [ok,msg] = checkSol(Opt);