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

加入变压器支路

This commit is contained in:
facat 2013-04-22 21:56:16 +08:00
parent 27bd5cb77c
commit 03865e7b64
9 changed files with 119 additions and 92 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
@SEOpti/SEOpti.m
View File

@ -20,6 +20,11 @@ classdef SEOpti
Balance=NaN;
mnle=NaN;
mnlrhs=NaN;
mTransP=NaN;
mTransQ=NaN;
transI=NaN;
transJ=NaN;
newwordParameter=NaN;
end
methods

2
@SEOpti/equ.m
View File

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

28
@SEOpti/fun.asv Normal file
View File

@ -0,0 +1,28 @@
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);%复数支路电流
SEBranchI=abs(cmpSEBranchI);% 支路电流幅值
Objective=Objective+sum((SEBranchI-this.mBranchI.^2).^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);
%% 变压器
newwordParameter=th
transP=TransPower( newwordParameter,cmpY,Volt,VAngel );
%% 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)));
end

22
@SEOpti/fun.m
View File

@ -1,6 +1,7 @@
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);%电压
@ -8,17 +9,24 @@ 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=real(cmpSEBranchI).^2+imag(cmpSEBranchI).^2;% ֧·µçÁ÷·ùÖµ
Objective=Objective+(SEBranchI-this.mBranchI.^2)'*(1./this.sigma^2*eye(length(this.mBranchI)))*(SEBranchI-this.mBranchI.^2);%%µçÁ÷
%% ֧·¹¦ÂÊ
SEBranchI=abs(cmpSEBranchI);% 支路电流幅值
Objective=Objective+sum((SEBranchI-this.mBranchI.^2).^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+(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);
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);
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.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+(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)));
end

6
@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)
function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ)
%INIT Summary of this function goes here
% Detailed explanation goes here
this.mVolt=mVolt;
@ -19,6 +19,10 @@ function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectio
this.mBranchI=mBranchI;
this.mBranchP=mBranchP;
this.mBranchQ=mBranchQ;
this.mTransP=mTransP;
this.mTransQ=mTransQ;
this.newwordParameter=newwordParameter;
% this.Y=Y;
% this.Yangle=Yangle;
% this.r=r;

12
TransPower.m Normal file
View File

@ -0,0 +1,12 @@
function [ output_args ] = TransPower( newwordParameter,cmpY,Volt,VAngel )
%TRANSPOWER Summary of this function goes here
% Detailed explanation goes here
transI=newwordParameter.trans.transI;
transJ=newwordParameter.trans.transJ;
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;
end

12
TransReactivePower.m Normal file
View File

@ -0,0 +1,12 @@
function [ output_args ] = TransReactivePower( newwordParameter,cmpY,Volt,VAngel )
%TRANSREAVTIVEPOWER Summary of this function goes here
% Detailed explanation goes here
transI=newwordParameter.trans.transI;
transJ=newwordParameter.trans.transJ;
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;
end

65
run.asv
View File

@ -1,8 +1,8 @@
clear
clc
yalmip('clear')
% yalmip('clear')
addpath('.\Powerflow')
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0');
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee14.dat', '0');
%% 量测量
% 电压 节点电流 支路电流 节点功率 支路功率
%%
@ -10,7 +10,7 @@ addpath('.\Powerflow')
% 电压 相角
%%
%% 开始生成量测量
sigma=0.05;% 标准差
sigma=1;% ±ê×¼²î
%% 电压
%电压幅值
rVolt=Volt; %幅值
@ -30,6 +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);%支路电流量测量
@ -45,7 +46,7 @@ rQD=QD;
QDi=find(QD~=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);
@ -53,61 +54,40 @@ 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(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG);
measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG);
fprintf('冗余度 %f\n',measurements/stateVarCount);
%% 进入状态估计计算
% SEVolt=sdpvar(length(Volt),1,'full');
% SEVAngel=sdpvar(length(Vangle),1,'full');
% Objective=(SEVolt-mVolt)'*(1./sigma^2*eye(length(mVolt)))*(SEVolt-mVolt);%电压
%% 支路电流
% cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压
% cmpSEBranchI=(cmpSEV(lineI)-cmpSEV(lineJ))./(lineR+1j*lineX);%复数支路电流
% SEBranchI=abs(cmpSEBranchI);% 支路电流幅值
% Objective=Objective+(SEBranchI-mBranchI)'*(1./sigma^2*eye(length(mBranchI)))*(SEBranchI-mBranchI);%%电流
%% 支路功率
% SEBranchP=real((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI));
% SEBranchQ=imag((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI));
% Objective=Objective+(SEBranchP-mBranchP)'*(1./sigma^2*eye(length(mBranchP)))*(SEBranchP-mBranchP);
% Objective=Objective+(SEBranchQ-mBranchQ)'*(1./sigma^2*eye(length(mBranchQ)))*(SEBranchQ-mBranchQ);
%% 0注入节点
zerosInjectionIndex=1:length(Volt);
zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
% PQ=diag(cmpSEV)*cmpY*cmpSEV;
% zeroInjP=real(PQ(zerosInjectionIndex));%% 0注入节点
% zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点
%% ±äѹÆ÷¹¦ÂÊ
rTransP=TransPower( newwordParameter,cmpY,rVolt,rVAngel );
rTransQ=TransReactivePower( newwordParameter,cmpY,rVolt,rVAngel );
mTransP=rTransP.*(normrnd(0,sigma,length(rTransP),1)+1);
mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1);
%% 发电机注入功率
% 先找到只有发电机的节点
PDQDi=union(PDi,QDi);
onlyPG=setdiff(PGi,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求解
% assign(SEVolt,1);
% assign(SEVAngel,1);
% Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0];
% Constraints=[SEVAngel(Balance)==0,PQ(zerosInjectionIndex)==0];
% options = sdpsettings('verbose',2,'solver','filtersd','usex0',1);
% solvesdp([],Objective,options)
% double(Objective)
% double(SEVolt)
% double(SEVAngel)
%% 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,mTransP,mTransQ);
opts = optiset('solver','ipopt');
opts.maxiter=2500;
opts.tolrfun=1e-3;
opts.tolafun=1e-3;
opts.warnings='warnings';
opts.display='iter';
opts.maxiter=85500;
opts.maxtime=3000;
opts.maxfeval=85000;
opts.maxnodes=85000;
opts.tolrfun=1e-4;
opts.tolafun=1e-4;
opts.warnings='all';
opts.display='off';
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)
% 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
%% 输出结果
@ -117,3 +97,4 @@ fprintf('Ŀ
fprintf('相对误差\n');
(abs(rVolt-double(SEVolt)))./(rVolt)
MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngel )

59
run.m
View File

@ -1,8 +1,8 @@
clear
clc
yalmip('clear')
% yalmip('clear')
addpath('.\Powerflow')
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0');
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee14.dat', '0');
%% 量测量
% 电压 节点电流 支路电流 节点功率 支路功率
%%
@ -10,7 +10,7 @@ addpath('.\Powerflow')
% 电压 相角
%%
%% 开始生成量测量
sigma=0.05;% 标准差
sigma=1;% ±ê×¼²î
%% 电压
%电压幅值
rVolt=Volt; %幅值
@ -30,6 +30,7 @@ lineJ=newwordParameter.line.lineJ;
lineR=newwordParameter.line.lineR;
lineX=newwordParameter.line.lineX;
lineB2=newwordParameter.line;
cmpBranchI=(cmpV(lineI)-cmpV(lineJ))./(lineR+1j*lineX);%复数支路电流
rBranchI=abs(cmpBranchI);% 支路电流幅值
mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量
@ -55,62 +56,38 @@ 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);
%% 进入状态估计计算
% SEVolt=sdpvar(length(Volt),1,'full');
% SEVAngel=sdpvar(length(Vangle),1,'full');
% Objective=(SEVolt-mVolt)'*(1./sigma^2*eye(length(mVolt)))*(SEVolt-mVolt);%电压
%% 支路电流
% cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压
% cmpSEBranchI=(cmpSEV(lineI)-cmpSEV(lineJ))./(lineR+1j*lineX);%复数支路电流
% SEBranchI=abs(cmpSEBranchI);% 支路电流幅值
% Objective=Objective+(SEBranchI-mBranchI)'*(1./sigma^2*eye(length(mBranchI)))*(SEBranchI-mBranchI);%%电流
%% 支路功率
% SEBranchP=real((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI));
% SEBranchQ=imag((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI));
% Objective=Objective+(SEBranchP-mBranchP)'*(1./sigma^2*eye(length(mBranchP)))*(SEBranchP-mBranchP);
% Objective=Objective+(SEBranchQ-mBranchQ)'*(1./sigma^2*eye(length(mBranchQ)))*(SEBranchQ-mBranchQ);
%% 0注入节点
zerosInjectionIndex=1:length(Volt);
zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
% PQ=diag(cmpSEV)*cmpY*cmpSEV;
% zeroInjP=real(PQ(zerosInjectionIndex));%% 0注入节点
% zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点
%% ±äѹÆ÷¹¦ÂÊ
rTransP=TransPower( newwordParameter,cmpY,rVolt,rVAngel );
rTransQ=TransReactivePower( newwordParameter,cmpY,rVolt,rVAngel );
mTransP=rTransP.*(normrnd(0,sigma,length(rTransP),1)+1);
mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1);
%% 发电机注入功率
% 先找到只有发电机的节点
PDQDi=union(PDi,QDi);
onlyPG=setdiff(PGi,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求解
% assign(SEVolt,1);
% assign(SEVAngel,1);
% Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0];
% Constraints=[SEVAngel(Balance)==0,PQ(zerosInjectionIndex)==0];
% options = sdpsettings('verbose',2,'solver','filtersd','usex0',1);
% solvesdp([],Objective,options)
% double(Objective)
% double(SEVolt)
% double(SEVAngel)
%% Opti ToolBox
Busnum=length(Volt);
seOpti=SEOpti();
seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ);
opts = optiset('solver','filtersd');
opts.maxiter=15500;
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;
opts.maxtime=3000;
opts.maxfeval=15000;
opts.maxnodes=15000;
opts.tolrfun=1e-5;
opts.tolafun=1e-5;
opts.maxfeval=85000;
opts.maxnodes=85000;
opts.tolrfun=1e-4;
opts.tolafun=1e-4;
opts.warnings='all';
opts.display='off';
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)
% 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
%% 输出结果