facat
/
stateestimation-ipm
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

用电流的平方而不是绝对值，因为abs的导数不好处理 

Signed-off-by: facat <dugg@21cn.com>
This commit is contained in:
facat 2013-04-14 22:17:01 +08:00
parent cf1cd0b95e
commit 059c5c2074
5 changed files with 13 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
@SEOpti/SEOpti.m
View File

@ -20,10 +20,6 @@ classdef SEOpti
Balance=NaN;
mnle=NaN;
mnlrhs=NaN;
Y=NaN;
Yangle=NaN;
r=NaN;
c=NaN;
end
methods

10
@SEOpti/equ.m
View File

@ -5,11 +5,11 @@ 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);
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);];
% 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);

4
@SEOpti/fun.m
View File

@ -8,8 +8,8 @@ Objective=(SEVolt-this.mVolt)'*(1./this.sigma^2*eye(length(this.mVolt)))*(SEVolt
% %% 支路电流
cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压
cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流
% SEBranchI=abs(cmpSEBranchI);% ֧·µçÁ÷·ùÖµ
% Objective=Objective+(SEBranchI-this.mBranchI)'*(1./this.sigma^2*eye(length(this.mBranchI)))*(SEBranchI-this.mBranchI);%%µçÁ÷
SEBranchI=real(cmpSEBranchI).^2+imag(cmpSEBranchI).^2;% ֧·µçÁ÷·ùÖµ
Objective=Objective+(SEBranchI-this.mBranchI.^2)'*(1./this.sigma^2*eye(length(this.mBranchI)))*(SEBranchI-this.mBranchI.^2);%%µçÁ÷
%% 支路功率
SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI));
SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI));

10
@SEOpti/init.m
View File

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

2
run.m
View File

@ -95,7 +95,7 @@ onlyQG=setdiff(QGi,PDQDi);
%% Opti ToolBox
Busnum=length(Volt);
seOpti=SEOpti();
seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,Y,Yangle,r,c );
seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ);
opts = optiset('solver','ipopt');
opts.maxiter=3500;
opts.tolrfun=1e-3;