加入了负荷电流

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

@Opti/Opti.asv

@@ -0,0 +1,34 @@
+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;
+        cl=NaN;%opti toolbox大于小于等于
+        cu=NaN;%等式、不等式约束右边的值
+        Balance=NaN;
+        rPD=NaN;
+        rQD=NaN;
+        m
+    end
+    
+    methods
+        function this=Opti()
+        end
+    end
+    
+end
+

@Opti/Opti.m

@@ -9,6 +9,7 @@ classdef Opti
         wPD=NaN;
         wQD=NaN;
         wVolt=NaN;
+        wLoadCurrent=NaN;
         PDi=NaN;
         QDi=NaN;
         Y=NaN;
@@ -22,6 +23,7 @@ classdef Opti
         Balance=NaN;
         rPD=NaN;
         rQD=NaN;
+        mLoadCurrent=NaN;
     end
     
     methods

@Opti/equ.m

@@ -23,7 +23,16 @@ 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
@@ -45,5 +54,6 @@ this.cl=[this.cl;0.8*rQD];
 output_args=full(output_args);
 this.cu=full(this.cu);
 this.cl=full(this.cl);
+
 end
 

@Opti/init.m

@@ -1,4 +1,4 @@
-function [ this ] = init( this,Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angel,r,c,PG,QG,Balance,rPD,rQD )
+function [ this ] = init( this,Volt0,PDi,QDi,wPD,wQD,wVolt,wLoadCurrent,PD0,QD0,rPD,rQD,Y,Angel,r,c,PG,QG,Balance,mLoadCurrent)
 %INIT Summary of this function goes here
 %   Detailed explanation goes here
 this.PDi=PDi;
@@ -17,5 +17,8 @@ this.Balance=Balance;
 this.Volt0=Volt0;
 this.rPD=rPD;
 this.rQD=rQD;
+this.wVolt=wVolt;
+this.wLoadCurrent=wLoadCurrent;
+this.mLoadCurrent=mLoadCurrent;
 end
 

@Opti/obj.asv

@@ -1,15 +1,22 @@
-function [ output_args ] = obj( this )
+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;
-PD=this.PD;
+% wVolt=this.wVolt;
-QD=this.QD;
+PD0=this.PD0;
-wVolt=t
+QD0=this.QD0;
-t4=wPD.*(PD(PD0~=0)-PD0(PD0~=0)).^2;
+t4=(wPD(PDi).*(PD-PD0(PDi))).^2;
-t5=wQD.*(QD(QD0~=0)-QD0(QD0~=0)).^2;
+t5=(wQD(QDi).*(QD-QD0(QDi))).^2;
-t6=wVolt.*(Volt-Volt0').^2;
+Volt0=this.Volt0;
-out_arg=sum(t4)+sum(t5)+sum(t6);
+SEVolt=x(length(PDi)+length(QDi)+1:length(PDi)+length(QDi)+length(Volt0));
-
+wVolt
+t6=(wVolt.*(SEVolt-Volt0)).^2;
+output_args=sum(t4)+sum(t5)+sum(t6);
+output_args=full(output_args);
 end
 

@Opti/obj.m

@@ -12,7 +12,19 @@ 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);
+Volt0=this.Volt0;
+SEVolt=x(length(PDi)+length(QDi)+1:length(PDi)+length(QDi)+length(Volt0));
+wVolt=this.wVolt;
+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;
+t7=(wLoadCurrent.*(loadCurrent-mLoadCurrent)).^2;
+output_args=sum(t4)+sum(t5)+sum(t7);
 output_args=full(output_args);
 end
 

LoadCurrent.asv

@@ -0,0 +1,6 @@
+function [ output_args ] = LoadCurrent( Volt,PD,QD,PDi,QDi )
+%% ¼ÆË㸺ºÉµçÁ÷
+
+
+end
+

LoadCurrent.m

@@ -0,0 +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)));
+output_args=t1;
+
+end
+

MaxDeviation.m

@@ -1,6 +1,8 @@
 function [output_arg]=MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD,rQD,PD,QD)
-t1=[rVolt;rVAngel;rPD;rQD];
+% t1=[rVolt;rVAngel;rPD;rQD];
-t2=[SEVolt;SEVAngel;PD;QD];
+% t2=[SEVolt;SEVAngel;PD;QD];
+t1=[rPD;rQD];
+t2=[PD;QD];
 t3=(t1(t1~=0)-double(t2(t1~=0)))./t1(t1~=0);
 t4=abs(t3);
 

OPF_Init.m

@@ -1,4 +1,4 @@
-function [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)
+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);
@@ -34,8 +34,11 @@ 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);
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
 % wD(11)=0;

Run_YALMIP.asv

@@ -1,11 +1,12 @@
 clc
 clear
-yalmip('clear')
+clear
+% yalmip('clear')
 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:/算例/东际911_2751267_2012-09-05/newFIle20.txt');
+    pf('E:/算例/新民Ⅰ906_2729823_2012-09-06/newFIle20.txt');
 
 %% ³±Á÷µÈʽ
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
@@ -18,7 +19,6 @@ PD0=PD;
 QD0=QD;
 PDReal=PD;%ÕæÖµ
 QDReal=QD;%ÕæÖµ
-%PD0(12)=PD0(12)+0.001;
 PG0(Balance)=PGBal(Balance);
 QG0(Balance)=QGBal(Balance);
 QG0(PVi)=QGBal(PVi);
@@ -36,19 +36,19 @@ xVolt=Volt;
 xUAngel=UAngel;
 % VMatrix=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
 % dP=PG-PD-diag(xVolt)*(Y.*cos(VMatrix))*xVolt';
-[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);
+rPD=PD;
-%% 定义变量
+rQD=QD;
+rVolt=Volt';
+rVAngel=UAngel';
 BalVolt=Volt(Balance);
-% Volt=sdpvar(Busnum,1);
+[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);
-% 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');
 %% ¼ÓÔëÒô
-load('PD0');
+load('20PD0');
-load('QD0');
+load('20QD0');
+mPD=PD0;
+mQD=QD0;
+mVolt=(1+normrnd(0,0.05,length(rVolt),1)).*rVolt;
 % PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0;
 % save('PD0','PD0');
 % QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0;
@@ -61,32 +61,46 @@ 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(Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle,r,c,PG,QG,Balance);
+seOpti=seOpti.init(mVolt,PDi,QDi,wPD,wQD,wVolt,wLoadCurrent,mPD,mQD,rPD(PDi),rQD(QDi),Y,Angle,r,c,PG,QG,Balance,mLoadCurrent);
-opts = optiset('solver','NLOPT');
+opts = optiset('solver','ipopt');
 opts.maxiter=85500;
-opts.maxtime=3000;
+opts.maxtime=30000;
 opts.maxfeval=85000;
 opts.maxnodes=85000;
 opts.tolrfun=1e-4;
 opts.tolafun=1e-4;
 opts.warnings='all';
 opts.display='iter';
-x0=[zeros(length(PDi)+length(QDi),1); ...
+x0=[0.8*rPD(PDi);0.8*rQD(QDi); ...
     ones(length(Volt),1); ...
     zeros(length(Volt),1)];
 % x0=[PD(PDi);QD(QDi);xVolt';xUAngel'];
 [~,seOpti]=seOpti.equ(x0);
-nlrhs=seOpti.Geteb();
+cl=seOpti.Getcl();
-nle=seOpti.Getgle();
+cu=seOpti.Getcu();
-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(Volt)*2+length(PDi)+length(QDi),'nl',@seOpti.equ,cl,cu,'options',opts)
-% Opt = opti('fun',@seOpti.obj,'ndec',length(PDi)+length(QDi),'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
 fval=seOpti.obj(x);
-fprintf('目标函数: %f\n',fval);
+fprintf('目标函数: %.20f\n',fval);
 toc
-rVolt=xVolt0';
+rVolt=Volt0';
 rVAngel=xUAngel';
-SEVolt=x(length(PDi)+length(QDi)+1:)
+SEVolt=x(length(PDi)+length(QDi)+1:length(Volt)+length(PDi)+length(QDi));
-MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
+SEVAngel=x(length(Volt)+length(PDi)+length(QDi)+1:end);
+fprintf('最大偏差\n')
+PD=x(1:length(PDi));
+QD=x(length(PDi)+1:length(PDi)+length(QDi));
+MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD(PDi),rQD(QDi),PD,QD)
+fprintf('统计偏差\n')
+StatDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD(PDi),rQD(QDi),PD,QD)
+%% 约束检查
+% seOpti.equ(x);
+% sum([SEVolt;PD;QD]>cu(length(SEVolt)*2+2:end));
+% sum([SEVolt;PD;QD]<cl(length(SEVolt)*2+2:end));

Run_YALMIP.m

@@ -1,11 +1,12 @@
 clc
 clear
-yalmip('clear')
+clear
+% yalmip('clear')
 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:/算例/东际911_2751267_2012-09-05/newFIle20.txt');
+    pf('E:/算例/新民Ⅰ906_2729823_2012-09-06/newFIle20.txt');
 
 %% ³±Á÷µÈʽ
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
@@ -18,7 +19,6 @@ PD0=PD;
 QD0=QD;
 PDReal=PD;%ÕæÖµ
 QDReal=QD;%ÕæÖµ
-%PD0(12)=PD0(12)+0.001;
 PG0(Balance)=PGBal(Balance);
 QG0(Balance)=QGBal(Balance);
 QG0(PVi)=QGBal(PVi);
@@ -38,19 +38,17 @@ xUAngel=UAngel;
 % dP=PG-PD-diag(xVolt)*(Y.*cos(VMatrix))*xVolt';
 rPD=PD;
 rQD=QD;
-[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);
+rVolt=Volt';
-%% 定义变量
+rVAngel=UAngel';
 BalVolt=Volt(Balance);
-% Volt=sdpvar(Busnum,1);
+[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);
-% 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');
 %% ¼ÓÔëÒô
-load('PD0');
+load('20PD0');
-load('QD0');
+load('20QD0');
+mPD=PD0;
+mQD=QD0;
+mVolt=(1+normrnd(0,0.05,length(rVolt),1)).*rVolt;
 % PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0;
 % save('PD0','PD0');
 % QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0;
@@ -63,30 +61,34 @@ 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(Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle,r,c,PG,QG,Balance,rPD(PDi),rQD(QDi));
+seOpti=seOpti.init(mVolt,PDi,QDi,wPD,wQD,wVolt,wLoadCurrent,mPD,mQD,rPD(PDi),rQD(QDi),Y,Angle,r,c,PG,QG,Balance,mLoadCurrent);
-opts = optiset('solver','filtersd');
+opts = optiset('solver','ipopt');
 opts.maxiter=85500;
-opts.maxtime=3000;
+opts.maxtime=30000;
 opts.maxfeval=85000;
 opts.maxnodes=85000;
 opts.tolrfun=1e-4;
 opts.tolafun=1e-4;
 opts.warnings='all';
-opts.display='iter';
+opts.display='off';
-% x0=[zeros(length(PDi)+length(QDi),1); ...
+x0=[0.8*rPD(PDi);0.8*rQD(QDi); ...
-%     ones(length(Volt),1); ...
+    ones(length(Volt),1); ...
-%     zeros(length(Volt),1)];
+    zeros(length(Volt),1)];
-x0=[PD(PDi);QD(QDi);xVolt';xUAngel'];
+% 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(PDi)+length(QDi),'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
 fval=seOpti.obj(x);
-fprintf('目标函数: %f\n',fval);
+fprintf('目标函数: %.20f\n',fval);
 toc
 rVolt=Volt0';
 rVAngel=xUAngel';
@@ -98,3 +100,7 @@ QD=x(length(PDi)+1:length(PDi)+length(QDi));
 MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD(PDi),rQD(QDi),PD,QD)
 fprintf('ͳ¼ÆÆ«²î\n')
 StatDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD(PDi),rQD(QDi),PD,QD)
+%% 约束检查
+% seOpti.equ(x);
+% sum([SEVolt;PD;QD]>cu(length(SEVolt)*2+2:end));
+% sum([SEVolt;PD;QD]<cl(length(SEVolt)*2+2:end));

StatDeviation.m

@@ -3,13 +3,14 @@ function [ output_args ] = StatDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD,rQD,P
 %   Detailed explanation goes here
 PDArray=PD;
 QDArray=QD;
+% QD0Array=rVAngel;
 PD0Array=repmat(rPD,1,1);
 QD0Array=repmat(rQD,1,1);
 PDDev=(PDArray-PD0Array)/0.05;
 QDDev=(QDArray-QD0Array)/0.05;
-% VoltDev=(SEVolt-rVolt)/0.05;
+VoltDev=(SEVolt-rVolt)/0.05;
-% VAngelDev=(SEVAngel-QD0Array)/0.05;
+VAngelDev=(SEVAngel-rVAngel)/0.05;
-wholeMat=[PDDev;QDDev;];
+wholeMat=[PDDev;QDDev;VoltDev;VAngelDev];
 t1=wholeMat.^2;
 t2=sum(t1,1);
 t3=t2/size(t1,1);