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

改了初值，效果很好. 

Signed-off-by: facat <dugg@21cn.com>
This commit is contained in:
facat 2013-04-23 11:43:35 +08:00
parent 03865e7b64
commit 202e7e62af
15 changed files with 81 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
@SEOpti/equ.m
View File

@ -10,9 +10,10 @@ PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV);
% 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.Balance+this.Busnum)];
out_arg=[out_arg;x(this.Balance+this.Busnum)];
out_arg=full(out_arg);
this.mnle=zeros(2*length(this.zerosInjectionIndex)+1,1);
this.mnle=zeros(length(out_arg),1);
this.mnlrhs=this.mnle;
end

15
@SEOpti/fun.asv
View File

@ -4,9 +4,8 @@ function [ Objective ] = fun(this, x )
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;
% %% ֧·µçÁ÷
Objective=sum(((SEVolt-this.mVolt)./this.sigma).^2);%µçѹ
%% ֧·µçÁ÷
cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压
cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流
SEBranchI=abs(cmpSEBranchI);% 支路电流幅值
@ -17,12 +16,16 @@ SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI));
Objective=Objective+sum((SEBranchP-this.mBranchP).^2/this.sigma^2);
Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/this.sigma^2);
%% 变压器
newwordParameter=th
transP=TransPower( newwordParameter,cmpY,Volt,VAngel );
newwordParameter=this.newwordParameter;
cmpY=this.cmpY;
transP=TransPower( newwordParameter,cmpY,SEVolt,SEVAngel );
Objective=Objective+sum((transP-this.mTransP).^2/this.sigma^2);
transQ=TransReactivePower( newwordParameter,cmpY,SEVolt,SEVAngel );
Objective=Objective+sum((transQ-this.mTransQ).^2/this.sigma^2);
%% 0注入节点
PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV);
%% 发电机注入功率
% 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.mPG(this.onlyPG)-real(PQ(this.onlyPG)))/this.sigmas.onlyPG)-real(PQ(this.onlyPG)));
% 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

21
@SEOpti/fun.m
View File

@ -4,29 +4,30 @@ function [ Objective ] = fun(this, x )
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;
% %% ֧·µçÁ÷
Objective=sum(((SEVolt-this.mVolt)./this.sigma).^2);%µçѹ
%% ֧·µçÁ÷
cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压
cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流
SEBranchI=abs(cmpSEBranchI);% 支路电流幅值
Objective=Objective+sum((SEBranchI-this.mBranchI.^2).^2./this.sigma^2);%%µçÁ÷
% (SEBranchI-this.mBranchI).^2./(this.sigma^2)
Objective=Objective+sum((SEBranchI-this.mBranchI).^2./(this.sigma^2));%%µçÁ÷
%% 线路功率
SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI));
SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI));
Objective=Objective+sum((SEBranchP-this.mBranchP).^2/this.sigma^2);
Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/this.sigma^2);
% (SEBranchP-this.mBranchP).^2/(this.sigma^2)./(this.mBranchP.^2);
Objective=Objective+sum((SEBranchP-this.mBranchP).^2/(this.sigma^2));
Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/(this.sigma^2));
%% 变压器
newwordParameter=this.newwordParameter;
cmpY=this.cmpY;
transP=TransPower( newwordParameter,cmpY,SEVolt,SEVAngel );
Objective=Objective+sum((transP-this.mTransP).^2/this.sigma^2);
Objective=Objective+sum((transP-this.mTransP).^2/(this.sigma^2));
transQ=TransReactivePower( newwordParameter,cmpY,SEVolt,SEVAngel );
Objective=Objective+sum((transQ-this.mTransQ).^2/this.sigma^2);
Objective=Objective+sum((transQ-this.mTransQ).^2/(this.sigma^2));
%% 0注入节点
PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV);
%% 发电机注入功率
% 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.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)));
% Objective=Objective+sum(((this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))/this.sigma).^2);
% Objective=Objective+sum(((this.mQG(this.onlyQG)-imag(PQ(this.onlyQG)))./this.sigma).^2);
end

2
TransPower.m
View File

@ -7,6 +7,6 @@ cmpSEV=Volt.*exp(1j*VAngel); %
% cmpY=sparse(r,c,Y.*(1j*YAngel),length(Volt),length(Volt));
cmpSETransI=(cmpSEV(transI)-cmpSEV(transJ)).*cmpY(sub2ind(size(cmpY),transI,transJ));
SETransP=real((cmpSEV(transI)-cmpSEV(transJ)).*conj(cmpSETransI));
output_args=SETransP;
output_args=full(SETransP);
end

2
TransReactivePower.m
View File

@ -7,6 +7,6 @@ cmpSEV=Volt.*exp(1j*VAngel); %
% cmpY=sparse(r,c,Y.*(1j*YAngel),length(Volt),length(Volt));
cmpSETransI=(cmpSEV(transI)-cmpSEV(transJ)).*cmpY(sub2ind(size(cmpY),transI,transJ));
SETransQ=imag((cmpSEV(transI)-cmpSEV(transJ)).*conj(cmpSETransI));
output_args=SETransQ;
output_args=full(SETransQ);
end

BIN
mBranchI.mat Normal file
View File

Binary file not shown.

BIN
mBranchP.mat Normal file
View File

Binary file not shown.

BIN
mBranchQ.mat Normal file
View File

Binary file not shown.

BIN
mPG.mat Normal file
View File

Binary file not shown.

BIN
mQG.mat Normal file
View File

Binary file not shown.

BIN
mTransP.mat Normal file
View File

Binary file not shown.

BIN
mTransQ.mat Normal file
View File

Binary file not shown.

BIN
mVolt.mat Normal file
View File

Binary file not shown.

38
run.asv
View File

@ -10,7 +10,7 @@ addpath('.\Powerflow')
% 电压 相角
%%
%% 开始生成量测量
sigma=1;% ±ê×¼²î
sigma=0.01;% ±ê×¼²î
%% 电压
%电压幅值
rVolt=Volt; %幅值
@ -30,7 +30,7 @@ lineJ=newwordParameter.line.lineJ;
lineR=newwordParameter.line.lineR;
lineX=newwordParameter.line.lineX;
lineB2=newwordParameter.line;
transI=
cmpBranchI=(cmpV(lineI)-cmpV(lineJ))./(lineR+1j*lineX);%复数支路电流
rBranchI=abs(cmpBranchI);% 支路电流幅值
mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量
@ -52,10 +52,6 @@ mPD=rPD.*(normrnd(0,sigma,length(rPD),1)+1);
mQD=rQD.*(normrnd(0,sigma,length(rQD),1)+1);
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(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG);
fprintf('ÈßÓà¶È %f\n',measurements/stateVarCount);
%% 0注入节点
zerosInjectionIndex=1:length(Volt);
zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
@ -69,9 +65,31 @@ mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1);
PDQDi=union(PDi,QDi);
onlyPG=setdiff(PGi,PDQDi);
onlyQG=setdiff(QGi,PDQDi);
%% ÈßÓà¶È¼ÆËã
stateVarCount=2*length(Volt);
measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG)+length(mTransP)+length(mTransQ);
fprintf('ÈßÓà¶È %f\n',measurements/stateVarCount);
%% Opti ToolBox
Busnum=length(Volt);
seOpti=SEOpti();
%% save
% save('mVolt','mVolt');
% save('mPG','mPG');
% save('mQG','mQG');
% save('mBranchI','mBranchI');
% save('mBranchP','mBranchP');
% save('mBranchQ','mBranchQ');
% save('mTransP','mTransP');
% save('mTransQ','mTransQ');
%% load
load('mVolt');
load('mPG');
load('mQG');
load('mBranchI');
load('mBranchP');
load('mBranchQ');
load('mTransP');
load('mTransQ');
seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ);
opts = optiset('solver','ipopt');
opts.maxiter=85500;
@ -82,12 +100,13 @@ opts.tolrfun=1e-4;
opts.tolafun=1e-4;
opts.warnings='all';
opts.display='off';
x0=[ones(Busnum,1);zeros(Busnum,1)];
% x0=[ones(Busnum,1);zeros(Busnum,1)];
x0=[rVolt;rVAngle];
[~,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)
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
%% 输出结果
@ -98,3 +117,4 @@ fprintf('
(abs(rVolt-double(SEVolt)))./(rVolt)
MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngel )
seOpti.fun(x);

36
run.m
View File

@ -10,7 +10,7 @@ addpath('.\Powerflow')
% 电压 相角
%%
%% 开始生成量测量
sigma=1;% ±ê×¼²î
sigma=0.03;% ±ê×¼²î
%% 电压
%电压幅值
rVolt=Volt; %幅值
@ -52,10 +52,6 @@ mPD=rPD.*(normrnd(0,sigma,length(rPD),1)+1);
mQD=rQD.*(normrnd(0,sigma,length(rQD),1)+1);
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(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG);
fprintf('ÈßÓàÈ %f\n',measurements/stateVarCount);
%% 0注入节点
zerosInjectionIndex=1:length(Volt);
zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
@ -69,9 +65,31 @@ mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1);
PDQDi=union(PDi,QDi);
onlyPG=setdiff(PGi,PDQDi);
onlyQG=setdiff(QGi,PDQDi);
%% ÈßÓàȼÆËã
stateVarCount=2*length(Volt);
measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG)+length(mTransP)+length(mTransQ);
fprintf('ÈßÓàÈ %f\n',measurements/stateVarCount);
%% Opti ToolBox
Busnum=length(Volt);
seOpti=SEOpti();
%% save
% save('mVolt','mVolt');
% save('mPG','mPG');
% save('mQG','mQG');
% save('mBranchI','mBranchI');
% save('mBranchP','mBranchP');
% save('mBranchQ','mBranchQ');
% save('mTransP','mTransP');
% save('mTransQ','mTransQ');
%% load
load('mVolt');
load('mPG');
load('mQG');
load('mBranchI');
load('mBranchP');
load('mBranchQ');
load('mTransP');
load('mTransQ');
seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ);
opts = optiset('solver','ipopt');
opts.maxiter=85500;
@ -82,12 +100,13 @@ opts.tolrfun=1e-4;
opts.tolafun=1e-4;
opts.warnings='all';
opts.display='off';
x0=[ones(Busnum,1);zeros(Busnum,1)];
x0=[ones(Busnum,1);-0.5*ones(Busnum,1)];
% x0=[rVolt;rVAngel];
[~,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)
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
%% 输出结果
@ -98,3 +117,4 @@ fprintf('
(abs(rVolt-double(SEVolt)))./(rVolt)
MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngel )
seOpti.fun(x);