不是用电流作为量测量，因为abs可导不好处理。

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

@SEOpti/SEOpti.m

@@ -1,24 +1,23 @@
 classdef SEOpti
     properties
-        mVolt=0;
+        mVolt=NaN;
-        mVAngel=0;
+        mBranchI=NaN;
-        mBranchI=0;
+        mBranchP=NaN;
-        mBranchP=0;
+        mBranchQ=NaN;
-        mBranchQ=0;
+        onlyPG=NaN;
-        onlyPG=0;
+        onlyQG=NaN;
-        onlyQG=0;
+        PG=NaN;
-        PG=0;
+        QG=NaN;
-        QG=0;
+        sigma=NaN;
-        sigma=0;
+        Busnum=NaN;
-        Busnum=0;
+        lineI=NaN;
-        lineI=0;
+        lineJ=NaN;
-        lineJ=0;
+        lineR=NaN;
-        lineR=0;
+        lineX=NaN;
-        lineX=0;
+        lineB2=NaN;
-        lineB2=0;
+        zerosInjectionIndex=NaN;
-        zerosInjectionIndex=0;
+        cmpY=NaN;
-        cmpY=0;
+        Balance=NaN;
-        Balance=0;
     end
     
     methods

@SEOpti/equ.m

@@ -5,7 +5,6 @@ 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);
-% Constraints=[SEVAngel(Balance)==0,PQ(zerosInjectionIndex)==0];
 out_arg=[real(PQ(this.zerosInjectionIndex));imag(PQ(this.zerosInjectionIndex));];
 
 out_arg=[out_arg;x(this.Busnum+this.Balance)];

@SEOpti/fun.m

@@ -1,23 +1,22 @@
 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);%电压
-    %% 支路电流
+% %% 支路电流
 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压
 cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流
-    SEBranchI=abs(cmpSEBranchI);% 支路电流幅值
+% SEBranchI=abs(cmpSEBranchI);% 支路电流幅值
-    Objective=Objective+(SEBranchI-this.mBranchI)'*(1./this.sigma^2*eye(length(this.mBranchI)))*(SEBranchI-this.mBranchI);%%电流
+% Objective=Objective+(SEBranchI-this.mBranchI)'*(1./this.sigma^2*eye(length(this.mBranchI)))*(SEBranchI-this.mBranchI);%%电流
 %% 支路功率
-    SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(SEBranchI));
+SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI));
-    SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(SEBranchI));
+SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI));
 Objective=Objective+(SEBranchP-this.mBranchP)'*(1./this.sigma^2*eye(length(this.mBranchP)))*(SEBranchP-this.mBranchP);
 Objective=Objective+(SEBranchQ-this.mBranchQ)'*(1./this.sigma^2*eye(length(this.mBranchQ)))*(SEBranchQ-this.mBranchQ);
-%     %% 0注入节点
+%% 0注入节点
 PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV);
-%     zeroInjP=real(PQ(this.zerosInjectionIndex));%% 0注入节点
-%     zeroInjQ=imag(PQ(this.zerosInjectionIndex));%% 0注入节点
 %% 发电机注入功率
 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.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)));

@SEOpti/init.asv

@@ -0,0 +1,24 @@
+function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance,mBranchI,mBranchP )
+%INIT Summary of this function goes here
+%   Detailed explanation goes here
+    this.mVolt=mVolt;
+    this.sigma=sigma;
+    this.Busnum=Busnum;
+    this.lineI=newwordParameter.line.lineI;
+    this.lineJ=newwordParameter.line.lineJ;
+    this.lineR=newwordParameter.line.lineR;
+    this.lineX=newwordParameter.line.lineX;
+    this.lineB2=newwordParameter.line;
+    this.zerosInjectionIndex=zerosInjectionIndex;
+    this.cmpY=cmpY;
+    this.onlyPG=onlyPG;
+    this.onlyQG=onlyQG;
+    this.PG=PG;
+    this.QG=QG;
+    this.Balance=Balance;
+    this.mBranchI=mBranchI;
+    this.mBranchP=mBranchP;
+    this.mBranchQ=mBranchQ;
+    
+end
+

@SEOpti/init.m

@@ -1,4 +1,4 @@
-function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance )
+function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance,mBranchI,mBranchP,mBranchQ )
 %INIT Summary of this function goes here
 %   Detailed explanation goes here
     this.mVolt=mVolt;
@@ -16,6 +16,9 @@ function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectio
     this.PG=PG;
     this.QG=QG;
     this.Balance=Balance;
+    this.mBranchI=mBranchI;
+    this.mBranchP=mBranchP;
+    this.mBranchQ=mBranchQ;
     
 end
 

run.m

@@ -10,7 +10,7 @@ addpath('.\Powerflow')
 % 电压 相角
 %%
 %% 开始生成量测量
-sigma=0.05;% ±ê×¼²î
+sigma=0.0005;% ±ê×¼²î
 %% 电压
 %电压幅值
 rVolt=Volt; %幅值
@@ -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');
@@ -95,12 +95,17 @@ 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);
+seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance,rBranchI,mBranchP,mBranchQ );
 nlrhs=seOpti.nlrhs();
 nle=seOpti.nle();
-Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle)
+opts = optiset('solver','ipopt');
-x0=ones(Busnum*2,1);
+opts.maxiter=2500;
+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)];
 [x,fval,exitflag,info] = solve(Opt,x0);
+info
 %% 输出结果
 SEVolt=x(1:length(Volt));
 SEVAngel=x(length(Volt)+Balance:2*length(Volt));