加了等式约束后很难收敛

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

@SEOpti/SEOpti.m

@@ -6,8 +6,8 @@ classdef SEOpti
         mBranchQ=NaN;
         onlyPG=NaN;
         onlyQG=NaN;
-        PG=NaN;
+        mPG=NaN;
-        QG=NaN;
+        mQG=NaN;
         sigma=NaN;
         Busnum=NaN;
         lineI=NaN;
@@ -18,6 +18,12 @@ classdef SEOpti
         zerosInjectionIndex=NaN;
         cmpY=NaN;
         Balance=NaN;
+        mnle=NaN;
+        mnlrhs=NaN;
+        Y=NaN;
+        Yangle=NaN;
+        r=NaN;
+        c=NaN;
     end
     
     methods

@SEOpti/equ.m

@@ -1,12 +1,18 @@
-function [ out_arg ] = equ(this, x )
+function [ out_arg,this ] = equ(this, x )
 %EQU Summary of this function goes here
 %   Detailed explanation goes here
 SEVolt=x(1:this.Busnum);
 SEVAngel=x(this.Busnum+1:2*this.Busnum);
 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压
 PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV);
-out_arg=[real(PQ(this.zerosInjectionIndex));imag(PQ(this.zerosInjectionIndex));];
+YAngle=sparse(this.r,this.c,SEVAngel(this.r)-SEVAngel(this.c)-this.Yangle,this.Busnum,this.Busnum);
-
+PQ1=diag(SEVolt)*(this.Y.*cos(YAngle))*SEVolt;
+PQ2=diag(SEVolt)*(this.Y.*sin(YAngle))*SEVolt;
+% out_arg=[real(PQ(this.zerosInjectionIndex));imag(PQ(this.zerosInjectionIndex));];
+out_arg=[PQ1(this.zerosInjectionIndex);PQ2(this.zerosInjectionIndex);];
 out_arg=[out_arg;x(this.Busnum+this.Balance)];
+out_arg=full(out_arg);
+this.mnle=zeros(2*length(this.zerosInjectionIndex)+1,1);
+this.mnlrhs=this.mnle;
 end
 

@SEOpti/fun.m

@@ -1,10 +1,10 @@
 function [ Objective ] = fun(this, x )
 %FUN Summary of this function goes here
 %   Detailed explanation goes here
-Objective=0;
 SEVolt=x(1:this.Busnum);
 SEVAngel=x(this.Busnum+1:2*this.Busnum);
 Objective=(SEVolt-this.mVolt)'*(1./this.sigma^2*eye(length(this.mVolt)))*(SEVolt-this.mVolt);%电压
+% Objective=0;
 % %% 支路电流
 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压
 cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流
@@ -18,7 +18,7 @@ Objective=Objective+(SEBranchQ-this.mBranchQ)'*(1./this.sigma^2*eye(length(this.
 %% 0注入节点
 PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV);
 %% 发电机注入功率
-Objective=Objective+(this.PG(this.onlyPG)-real(PQ(this.onlyPG)))'*(1./this.sigma^2*eye(length(this.PG(this.onlyPG))))*(this.PG(this.onlyPG)-real(PQ(this.onlyPG)));
+Objective=Objective+(this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))'*(1./this.sigma^2*eye(length(this.mPG(this.onlyPG))))*(this.mPG(this.onlyPG)-real(PQ(this.onlyPG)));
-Objective=Objective+(this.QG(this.onlyQG)-imag(PQ(this.onlyQG)))'*(1./this.sigma^2*eye(length (this.QG(this.onlyQG))))*(this.QG(this.onlyQG)-imag(PQ(this.onlyQG)));
+Objective=Objective+(this.mQG(this.onlyQG)-imag(PQ(this.onlyQG)))'*(1./this.sigma^2*eye(length(this.mQG(this.onlyQG))))*(this.mQG(this.onlyQG)-imag(PQ(this.onlyQG)));
 end
 

@SEOpti/init.m

@@ -1,4 +1,4 @@
-function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance,mBranchI,mBranchP,mBranchQ )
+function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,Y,Yangle,r,c )
 %INIT Summary of this function goes here
 %   Detailed explanation goes here
     this.mVolt=mVolt;
@@ -13,12 +13,15 @@ function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectio
     this.cmpY=cmpY;
     this.onlyPG=onlyPG;
     this.onlyQG=onlyQG;
-    this.PG=PG;
+    this.mPG=mPG;
-    this.QG=QG;
+    this.mQG=mQG;
     this.Balance=Balance;
     this.mBranchI=mBranchI;
     this.mBranchP=mBranchP;
     this.mBranchQ=mBranchQ;
-    
+    this.Y=Y;
+    this.Yangle=Yangle;
+    this.r=r;
+    this.c=c;
 end
 

@SEOpti/nle.m

@@ -1,7 +1,7 @@
 function [ output_args ] = nle( this )
 %NLE Summary of this function goes here
 %   Detailed explanation goes here
-output_args=zeros(2*length(this.zerosInjectionIndex)+1,1);
+output_args=this.mnle;
 
 end
 

@SEOpti/nlrhs.m

@@ -1,7 +1,7 @@
 function [ output_args ] = nlrhs( this )
 %NLRHS Summary of this function goes here
 %   Detailed explanation goes here
-output_args=zeros(2*length(this.zerosInjectionIndex)+1,1);
+output_args=this.mnlrhs;
 
 end
 

run.asv

@@ -2,7 +2,7 @@ clear
 clc
 yalmip('clear')
 addpath('.\Powerflow')
-[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee4.dat', '0');
+[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0');
 %% 量测量
 % 电压 节点电流 支路电流 节点功率 支路功率
 %%
@@ -39,11 +39,11 @@ mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%֧·
 rBranchQ=imag((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI));
 mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量
 %% 注入功率
-rPD=PD(PD~=0);
+rPD=PD;
 PDi=find(PD~=0);
-rQD=QD(QD~=0);
+rQD=QD;
 QDi=find(QD~=0);
-rPG=PG(PG~=0);
+rPG=PG;
 PGi=find(PG~=0);
 rQG=QG(QG~=0);
 QGi=find(QG~=0);
@@ -53,7 +53,7 @@ mPG=rPG.*(normrnd(0,sigma,length(rPG),1)+1);
 mQG=rQG.*(normrnd(0,sigma,length(rQG),1)+1);
 %% 冗余度计算
 stateVarCount=2*length(Volt);
-measurements=length(mVolt)+length(mI)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG);
+measurements=length(mVolt)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG);
 fprintf('冗余度 %f\n',measurements/stateVarCount);
 %% 进入状态估计计算
 % SEVolt=sdpvar(length(Volt),1,'full');
@@ -77,9 +77,9 @@ zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
 % zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点
 %% 发电机注入功率
 % 先找到只有发电机的节点
-% PDQDi=union(PDi,QDi);
+PDQDi=union(PDi,QDi);
-% onlyPG=setdiff(PGi,PDQDi);
+onlyPG=setdiff(PGi,PDQDi);
-% onlyQG=setdiff(QGi,PDQDi);
+onlyQG=setdiff(QGi,PDQDi);
 % Objective=Objective+(PG(onlyPG)-real(PQ(onlyPG)))'*(1./sigma^2*eye(length(PG(onlyPG))))*(PG(onlyPG)-real(PQ(onlyPG)));
 % Objective=Objective+(QG(onlyQG)-imag(PQ(onlyQG)))'*(1./sigma^2*eye(length(QG(onlyQG))))*(QG(onlyQG)-imag(PQ(onlyQG)));
 %% YALMIP求解
@@ -95,13 +95,25 @@ zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
 %% Opti ToolBox
 Busnum=length(Volt);
 seOpti=SEOpti();
-seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex );
+seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,Y,Yangle,r,c );
-Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2)
+opts = optiset('solver','ipopt');
-x0=ones(Busnum*2,1);
+opts.maxiter=2500;
+opts.tolrfun=1e-3;
+opts.tolafun=1e-3;
+opts.warnings='warnings';
+opts.display='iter';
+x0=[ones(Busnum,1);zeros(Busnum,1)];
+[~,seOpti]=seOpti.equ(x0);
+nlrhs=seOpti.nlrhs();
+nle=seOpti.nle();
+Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts)
+% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts)
 [x,fval,exitflag,info] = solve(Opt,x0);
 info
-x
 %% 输出结果
-% fprintf('Ïà¶ÔÎó²î\n');
+SEVolt=x(1:length(Volt));
-% (abs(rVolt-double(SEVolt)))./(rVolt)
+SEVAngel=x(length(Volt)+Balance:2*length(Volt));
-% MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
+fprintf('Ä¿±êº¯ÊýΪ:%f\n',fval);
+fprintf('Ïà¶ÔÎó²î\n');
+(abs(rVolt-double(SEVolt)))./(rVolt)
+MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)

run.m

@@ -2,7 +2,7 @@ clear
 clc
 yalmip('clear')
 addpath('.\Powerflow')
-[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee14.dat', '0');
+[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0');
 %% 量测量
 % 电压 节点电流 支路电流 节点功率 支路功率
 %%
@@ -10,7 +10,7 @@ addpath('.\Powerflow')
 % 电压 相角
 %%
 %% 开始生成量测量
-sigma=0.0005;% ±ê×¼²î
+sigma=0.05;% ±ê×¼²î
 %% 电压
 %电压幅值
 rVolt=Volt; %幅值
@@ -39,13 +39,13 @@ mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%֧·
 rBranchQ=imag((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI));
 mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量
 %% 注入功率
-rPD=PD(PD~=0);
+rPD=PD;
 PDi=find(PD~=0);
-rQD=QD(QD~=0);
+rQD=QD;
 QDi=find(QD~=0);
-rPG=PG(PG~=0);
+rPG=PG;
 PGi=find(PG~=0);
-rQG=QG(QG~=0);
+rQG=QG;
 QGi=find(QG~=0);
 mPD=rPD.*(normrnd(0,sigma,length(rPD),1)+1);
 mQD=rQD.*(normrnd(0,sigma,length(rQD),1)+1);
@@ -95,15 +95,19 @@ onlyQG=setdiff(QGi,PDQDi);
 %% Opti ToolBox
 Busnum=length(Volt);
 seOpti=SEOpti();
-seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance,rBranchI,mBranchP,mBranchQ );
+seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,Y,Yangle,r,c );
-nlrhs=seOpti.nlrhs();
-nle=seOpti.nle();
 opts = optiset('solver','ipopt');
 opts.maxiter=2500;
+opts.tolrfun=1e-3;
+opts.tolafun=1e-3;
+opts.warnings='warnings';
 opts.display='iter';
-% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts)
-Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts)
 x0=[ones(Busnum,1);zeros(Busnum,1)];
+[~,seOpti]=seOpti.equ(x0);
+nlrhs=seOpti.nlrhs();
+nle=seOpti.nle();
+Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts)
+% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts)
 [x,fval,exitflag,info] = solve(Opt,x0);
 info
 %% 输出结果