改成了含0注入约束（不是潮流约束）的最小二乘法。

Signed-off-by: dmy@lab <dmy@lab.lab>

@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));
+Volt=x(1:length(Volt0));
-VAngel=x(length(PDi)+length(QDi)+length(Volt0)+1:end);
+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);

@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);
+VMatrix=sparse(r,c,SEVAngel(r)-SEVAngel(c)-Angel,busNum,busNum);
-dP=PG-PD_-diag(Volt)*(Y.*cos(VMatrix))*Volt;
+zerosP=diag(SEVolt)*(Y.*cos(VMatrix))*SEVolt;
-dQ=QG-QD_-diag(Volt)*(Y.*spfun(@sin,VMatrix))*Volt;
+zerosQ=diag(SEVolt)*(Y.*spfun(@sin,VMatrix))*SEVolt;
-output_args=[dP;dQ];%潮流约束
+PDi=this.PDi;
-%% 使用零注入begin
+QDi=this.QDi;
-% zeroInj=setdiff(1:busNum,union(Balance,union(PDi,QDi)));
+zeroInj=setdiff(1:busNum,union(Balance,union(PDi,QDi)));
-% output_args=[output_args(zeroInj);output_args(length(Volt0)+zeroInj)];
+output_args=[zerosP(zeroInj);zerosQ(zeroInj)];
-%% 使用零注入end
+this.cu=zeros(length(output_args),1);
-%电压恒定
+this.cl=zeros(length(output_args),1);
-output_args=[output_args;Volt(Balance)-1];
+%% 使用零注入end
-%相角恒定
+%电压恒定
-output_args=[output_args;VAngel(Balance)];
+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);

@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))).^2;
-% 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);
+SEVAngel=x(length(SEVolt)+1:end);
-loadCurrent=LoadCurrent( SEVolt,SEVAngel,PD0(PDi),QD0(QDi),PDi,QDi );
+busNum=length(Volt0);
-wLoadCurrent=this.wLoadCurrent;
+Angel=this.Angel;
-mLoadCurrent=this.mLoadCurrent;
+r=this.r;
-t7=(wLoadCurrent.*(loadCurrent-mLoadCurrent)).^2;
+c=this.c;
-t7=wLoadCurrent.*((loadCurrent-mLoadCurrent)./mLoadCurrent).^2;
+PG=this.PG;
-output_args=sum(t4)+sum(t5)+sum(t6)+sum(0);
+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
 

@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);
+SEVAngel=x(length(SEVolt)+1:end);
-loadCurrent=LoadCurrent( SEVolt,SEVAngel,PD0(PDi),QD0(QDi),PDi,QDi );
+busNum=length(Volt0);
-wLoadCurrent=this.wLoadCurrent;
+Angel=this.Angel;
-mLoadCurrent=this.mLoadCurrent;
+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;
+% t7=wLoadCurrent.*((loadCurrent-mLoadCurrent)./mLoadCurrent).^2;
-output_args=sum(t4)+sum(t5)+sum(t6);
+% output_args=sum(t4)+sum(t5)+sum(t6);
+output_args=sum(t6)+sum(t4)+sum(t5);
 output_args=full(output_args);
+
+
+
 end
 

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;
+% wPD(1:end)=0;
-wQD(1:end)=0;
+% wQD(1:end)=0;
-wVolt(1:end)=0;
+% wVolt(1:end)=0;
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
 % wD(11)=0;

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); ...
+    x0=[ones(length(Volt),1); ...
-        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));
+    SEVolt=x(1:length(Volt));
-    SEVAngel=x(length(Volt)+length(PDi)+length(QDi)+1:end);
+    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);