diff --git a/@Opti/equ.asv b/@Opti/equ.asv index 2c640b2..d902550 100644 --- a/@Opti/equ.asv +++ b/@Opti/equ.asv @@ -23,15 +23,37 @@ 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;VAngel(Balance)]; +output_args=[dP;dQ]; +%% 使用零注入begin +% zeroInj=setdiff(1:busNum,union(Balance,union(PDi,QDi))); +% output_args=[output_args(zeroInj);output_args(length(PDi)+zeroInj)]; +%% 使用零注入end + +%电压恒定 +output_args=[output_args;Volt(Balance)-1]; +%相角恒定 +output_args=[output_args;VAngel(Balance)]; +this.cu=zeros(length(output_args),1); +this.cl=this.cu; %% 开始增加不等式约束-电压\PD\QD %电压不等式约束 output_args=[output_args;Volt]; +this.cu=[this.cu;1.1*ones(length(Volt),1)]; +this.cl=[this.cl;0.9*ones(length(Volt),1)]; +%% PD +rPD=this.rPD; +output_args=[output_args;PD]; +this.cu=[this.cu;1.5*rPD]; +this.cl=[this.cl;0.5*rPD]; +%% QD +rQD=this.rQD; +output_args=[output_args;QD]; +this.cu=[this.cu;1.9*rQD]; +this.cl=[this.cl;0.9*rQD]; +%% 稠密化 output_args=full(output_args); -this.cu=zeros(length(output_args),1); -this.cl=this.cu; -%电压 -this.cu=[this.cu;1.1]; -this.cl=this.cu; +this.cu=full(this.cu); +this.cl=full(this.cl); + end diff --git a/@Opti/equ.m b/@Opti/equ.m index 0ad51e0..6ef94c3 100644 --- a/@Opti/equ.m +++ b/@Opti/equ.m @@ -26,7 +26,7 @@ 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(PDi)+zeroInj)]; +% output_args=[output_args(zeroInj);output_args(length(Volt0)+zeroInj)]; %% 使用零注入end %电压恒定 @@ -38,8 +38,8 @@ this.cl=this.cu; %% 开始增加不等式约束-电压\PD\QD %电压不等式约束 output_args=[output_args;Volt]; -this.cu=[this.cu;1.1*ones(length(Volt),1)]; -this.cl=[this.cl;0.9*ones(length(Volt),1)]; +this.cu=[this.cu;1.07*ones(length(Volt),1)]; +this.cl=[this.cl;0.93*ones(length(Volt),1)]; %% PD rPD=this.rPD; output_args=[output_args;PD]; diff --git a/@Opti/obj.asv b/@Opti/obj.asv index 5a5d070..24ca46b 100644 --- a/@Opti/obj.asv +++ b/@Opti/obj.asv @@ -7,16 +7,28 @@ 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; +% 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)); -wVolt +wVolt=this.wVolt; t6=(wVolt.*(SEVolt-Volt0)).^2; -output_args=sum(t4)+sum(t5)+sum(t6); +% t6=wVolt.*((SEVolt-Volt0)./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); output_args=full(output_args); end diff --git a/@Opti/obj.m b/@Opti/obj.m index 8597e41..9fb5d67 100644 --- a/@Opti/obj.m +++ b/@Opti/obj.m @@ -7,24 +7,30 @@ 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; +% 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)); wVolt=this.wVolt; -t6=(wVolt.*(SEVolt-Volt0)).^2; +% t6=(wVolt.*(SEVolt-Volt0)).^2; +t6=wVolt.*((SEVolt-Volt0)./Volt0).^2; Balance=this.Balance; t6(Balance)=0;%不计算平衡节点的 +% xx=t6(8); +% t6=xx; +% t6(9:end)=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; -output_args=sum(t4)+sum(t5)+sum(t7); +% t7=(wLoadCurrent.*(loadCurrent-mLoadCurrent)).^2; +t7=wLoadCurrent.*((loadCurrent-mLoadCurrent)./mLoadCurrent).^2; +output_args=sum(t4)+sum(t5)+sum(t6)+sum(0); output_args=full(output_args); end diff --git a/LoadCurrent.asv b/LoadCurrent.asv index 924bfc1..cbfaffc 100644 --- a/LoadCurrent.asv +++ b/LoadCurrent.asv @@ -1,6 +1,9 @@ -function [ output_args ] = LoadCurrent( Volt,PD,QD,PDi,QDi ) +function [ output_args ] = LoadCurrent( Volt,VAngel,PD,QD,PDi,QDi ) %% 计算负荷电流 - +cmpVolt=Volt.*exp(1j*VAngel); +% t1=real(conj((PD+1j*QD)./cmpVolt(PDi))); +t1=real(conj((PD+1j*QD))./); +output_args=t1; end diff --git a/LoadCurrent.m b/LoadCurrent.m index 7f74a65..2d617c6 100644 --- a/LoadCurrent.m +++ b/LoadCurrent.m @@ -1,7 +1,8 @@ function [ output_args ] = LoadCurrent( Volt,VAngel,PD,QD,PDi,QDi ) %% 计算负荷电流 cmpVolt=Volt.*exp(1j*VAngel); -t1=abs(conj((PD+1j*QD)./cmpVolt(PDi))); +t1=real(conj((PD+1j*QD)./cmpVolt(PDi))); +% t1=real(conj((PD+1j*QD)).*Volt(PDi)); output_args=t1; end diff --git a/MaxDeviation.m b/MaxDeviation.m index 3c7a897..6cd0287 100644 --- a/MaxDeviation.m +++ b/MaxDeviation.m @@ -5,7 +5,6 @@ t1=[rPD;rQD]; t2=[PD;QD]; t3=(t1(t1~=0)-double(t2(t1~=0)))./t1(t1~=0); t4=abs(t3); - output_arg=max(t4); end diff --git a/OPF_Init.asv b/OPF_Init.asv new file mode 100644 index 0000000..375b0f2 --- /dev/null +++ b/OPF_Init.asv @@ -0,0 +1,52 @@ +function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,wLoadCurrent,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD) +Loadi=find(QD~=0 | PD~=0); +PDi=find(PD~=0); +QDi=find(QD~=0); +notLoadi=setdiff(1:Busnum,Loadi); +PGOnly=setdiff(PGi,PDi); +QGOnly=setdiff(PVi,QDi); +PDOnly=setdiff(PDi,PGi); +QDOnly=setdiff(QDi,PVi); +PGPD=intersect(PDi,PGi); +QGQD=intersect(QDi,PVi); +noPGPDQGQD=setdiff(1:Busnum,union(union(PDi,PGi),union(QDi,PVi))); +%Loadi=[1:Busnum]'; +RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*1+Busnum*1; %约束条件数,放开所有QD,不单独约束QD上下限 +RestraintCount=RestraintCount+length(Loadi); %加上功率因数,用功率因数约束QD +t_Bal_volt=Volt(Balance); +Volt=sparse(1*ones(1,Busnum)); +Volt(Balance)=t_Bal_volt; +UAngel=sparse(1,Busnum); +Init_Z=sparse(ones(1,RestraintCount)); +Init_W=sparse(-1*ones(1,RestraintCount)); +Init_L=1*sparse(ones(1,RestraintCount)); +Init_U=1*sparse(ones(1,RestraintCount)); +Init_Y=sparse(1,2*Busnum);%与学姐一致 +tPU=sparse(GenU(:,2));% 发电机有功上限 +tQU=sparse(PVQU(:,1));% 无功上限 +tPL=sparse(GenL(:,2));% 发电机有功下限 +tQL=sparse(PVQL(:,1));% 无功下限 +% PG(PGi)=(tPU+tPL)/2; +% QG(PVi)=(tQU+tQL)/2; +wPG=1*ones(size(PGi,1),1); +wQG=1*ones(size(PVi,1),1); +%randInt=randperm(size(Loadi,1)); +%randPDind=randInt(1:10); +randPDind=0; +wPD=1/0.05*zeros(length(PD),1); +wPD(1:2:end)=0; +wQD=1/0.05*zeros(length(QD),1); +wQD(1:2:end)=0; +wVolt=1/0.05*zeros(Busnum,1); +wVolt(1:2:end)=0; +wLoadCurrent=1/0.05*zeros(length(PDi),1); +% wLoadCurrent(1:2:end)=0; +%wD(randPDind)=0;%一些负荷不约束 +%wD(7)=0; +% wD(11)=0; +PD=1*PD0; +%powerFacter=0.98; +%QD=PD*sqrt((1-powerFacter^2)/powerFacter^2); +QD=QD0; + +end \ No newline at end of file diff --git a/OPF_Init.m b/OPF_Init.m index 4e66492..beaf250 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -34,11 +34,24 @@ wQG=1*ones(size(PVi,1),1); %randPDind=randInt(1:10); randPDind=0; wPD=1/0.05*ones(length(PD),1); -% wPD(1:2:end)=0; wQD=1/0.05*ones(length(QD),1); -% wQD(1:2:end)=0; wVolt=1/0.05*ones(Busnum,1); wLoadCurrent=1/0.05*ones(length(PDi),1); +%% 70% +% noLoad=[1;4;5;7;8;10;11;13;14;16;17;19;20;21]; +% wPD(noLoad)=0; +% wQD(noLoad)=0; +% wVolt(noLoad)=0; +%% 50% +% wPD(1:2:end)=0; +% wQD(1:2:end)=0; +% wVolt(1:2:end)=0; +% wLoadCurrent(1:2:end)=0; +%% 30% +noLoad=[1;5;7;11;15;20]; +wPD(noLoad)=0; +wQD(noLoad)=0; +wVolt(noLoad)=0; %wD(randPDind)=0;%一些负荷不约束 %wD(7)=0; % wD(11)=0; diff --git a/PD0.mat b/PD0.mat index 39a84c6..2b451de 100644 Binary files a/PD0.mat and b/PD0.mat differ diff --git a/QD0.mat b/QD0.mat index ae05a93..6a15376 100644 Binary files a/QD0.mat and b/QD0.mat differ diff --git a/Run_YALMIP.m b/Run_YALMIP.m index 6ebc6df..4df1d06 100644 --- a/Run_YALMIP.m +++ b/Run_YALMIP.m @@ -6,7 +6,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('E:/算例/新民Ⅰ906_2729823_2012-09-06/newFIle20.txt'); + pf('E:/算例/东际911_2751267_2012-09-05/newFIle20.txt'); %% 潮流等式 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); @@ -42,31 +42,37 @@ rVolt=Volt'; rVAngel=UAngel'; BalVolt=Volt(Balance); [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,wLoadCurrent,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD); - +%% 开始 +Busnum=length(Volt); +PDi=find(PD~=0); +QDi=find(QD~=0); +%% 电流真实值 +rLoadCurrent=LoadCurrent( rVolt,rVAngel,rPD(PDi),rQD(QDi),PDi,QDi ); %% 加噪音 -load('20PD0'); -load('20QD0'); -mPD=PD0; -mQD=QD0; -mVolt=(1+normrnd(0,0.05,length(rVolt),1)).*rVolt; +load('PD0'); +load('QD0'); +load('mVolt'); +load('mLoadCurrent'); +% mVolt=(1+normrnd(0,0.05,length(rVolt),1)).*rVolt; +% save('mVolt','mVolt'); % PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0; % save('PD0','PD0'); % QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0; % save('QD0','QD0'); +%% 电流测量值 +% mLoadCurrent=(1+normrnd(0,0.05,length(rLoadCurrent),1)).*rLoadCurrent; +% save('mLoadCurrent','mLoadCurrent'); +mPD=PD0; +mQD=QD0; %% 目标函数 %% Opti Toolbox -Busnum=length(Volt); -PDi=find(PD~=0); -QDi=find(QD~=0); + % PD0=PD0(PDi); % QD0=QD0(QDi); -%% 电流真实值 -rLoadCurrent=LoadCurrent( rVolt,rVAngel,rPD(PDi),rQD(QDi),PDi,QDi ); -%% 电流测量值 -mLoadCurrent=(1+normrnd(0,0.05,length(rLoadCurrent),1)).*rLoadCurrent; + seOpti=Opti(); -seOpti=seOpti.init(mVolt,PDi,QDi,wPD,wQD,wVolt,wLoadCurrent,mPD,mQD,rPD(PDi),rQD(QDi),Y,Angle,r,c,PG,QG,Balance,mLoadCurrent); +seOpti=seOpti.init(rVolt*1.005,PDi,QDi,wPD,wQD,wVolt,wLoadCurrent,mPD,mQD,rPD(PDi),rQD(QDi),Y,Angle,r,c,PG,QG,Balance,mLoadCurrent); opts = optiset('solver','ipopt'); opts.maxiter=85500; opts.maxtime=30000; @@ -76,7 +82,7 @@ opts.tolrfun=1e-4; opts.tolafun=1e-4; opts.warnings='all'; opts.display='off'; -x0=[0.8*rPD(PDi);0.8*rQD(QDi); ... +x0=[0.99*rPD(PDi);0.99*rQD(QDi); ... ones(length(Volt),1); ... zeros(length(Volt),1)]; % x0=[PD(PDi);QD(QDi);xVolt';xUAngel']; diff --git a/StatDeviation.m b/StatDeviation.m index 0818987..1af90ff 100644 --- a/StatDeviation.m +++ b/StatDeviation.m @@ -10,7 +10,7 @@ PDDev=(PDArray-PD0Array)/0.05; QDDev=(QDArray-QD0Array)/0.05; VoltDev=(SEVolt-rVolt)/0.05; VAngelDev=(SEVAngel-rVAngel)/0.05; -wholeMat=[PDDev;QDDev;VoltDev;VAngelDev]; +wholeMat=[VoltDev;VAngelDev]; t1=wholeMat.^2; t2=sum(t1,1); t3=t2/size(t1,1); diff --git a/mLoadCurrent.mat b/mLoadCurrent.mat new file mode 100644 index 0000000..602824d Binary files /dev/null and b/mLoadCurrent.mat differ diff --git a/mVolt.mat b/mVolt.mat new file mode 100644 index 0000000..b2562f5 Binary files /dev/null and b/mVolt.mat differ