改了目标函数的形式

Signed-off-by: facat <dugg@21cn.com>

@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.cu=full(this.cu);
-this.cl=this.cu;
+this.cl=full(this.cl);
-%µçѹ
+
-this.cu=[this.cu;1.1];
-this.cl=this.cu;
 end
 

@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.cu=[this.cu;1.07*ones(length(Volt),1)];
-this.cl=[this.cl;0.9*ones(length(Volt),1)];
+this.cl=[this.cl;0.93*ones(length(Volt),1)];
 %% PD
 rPD=this.rPD;
 output_args=[output_args;PD];

@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
 

@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;
+% t4=(wPD(PDi).*(PD-PD0(PDi))).^2;
-t5=(wQD(QDi).*(QD-QD0(QDi))).^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;
+% t7=(wLoadCurrent.*(loadCurrent-mLoadCurrent)).^2;
-output_args=sum(t4)+sum(t5)+sum(t7);
+t7=wLoadCurrent.*((loadCurrent-mLoadCurrent)./mLoadCurrent).^2;
+output_args=sum(t4)+sum(t5)+sum(t6)+sum(0);
 output_args=full(output_args);
 end
 

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
 

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

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
 

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

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;

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('PD0');
-load('20QD0');
+load('QD0');
-mPD=PD0;
+load('mVolt');
-mQD=QD0;
+load('mLoadCurrent');
-mVolt=(1+normrnd(0,0.05,length(rVolt),1)).*rVolt;
+% 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'];

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);