stateestimation-self-derivatives-measurementdistribution/run.m

clear
clc
yalmip('clear')
addpath('.\Powerflow')
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0');
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD>ѹ <EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> ֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD> <EFBFBD>ڵ㹦<EFBFBD><EFBFBD> ֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%%
%% ״̬<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD>ѹ <EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%%
%% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
sigma=0.05;% <EFBFBD><EFBFBD>׼<EFBFBD><EFBFBD>
%% <EFBFBD><EFBFBD>ѹ
%<EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD>ֵ
rVolt=Volt; %<EFBFBD><EFBFBD>ֵ
mVolt=rVolt.*(normrnd(0,sigma,length(Volt),1)+1);%<EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
rVAngel=Vangle;
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%ע<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
cmpY=Y.*exp(1j*sparse(r,c,Yangle,length(Y),length(Y)));%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɾ<EFBFBD><EFBFBD><EFBFBD>
cmpV=Volt.*exp(1j*Vangle); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
cmpI=cmpY*cmpV;% ע<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
rI=abs(cmpI); %ע<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҫ<EFBFBD><EFBFBD><EFBFBD>ǵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
mI=rI.*(normrnd(0,sigma,length(rI),1)+1);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% ֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% ֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
lineI=newwordParameter.line.lineI;
lineJ=newwordParameter.line.lineJ;
lineR=newwordParameter.line.lineR;
lineX=newwordParameter.line.lineX;
lineB2=newwordParameter.line;
cmpBranchI=(cmpV(lineI)-cmpV(lineJ))./(lineR+1j*lineX);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
rBranchI=abs(cmpBranchI);% ֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% ֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
rBranchP=real((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI));
mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
rBranchQ=imag((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI));
mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% ע<EFBFBD>빦<EFBFBD><EFBFBD>
rPD=PD;
PDi=find(PD~=0);
rQD=QD;
QDi=find(QD~=0);
rPG=PG;
PGi=find(PG~=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);
mPG=rPG.*(normrnd(0,sigma,length(rPG),1)+1);
mQG=rQG.*(normrnd(0,sigma,length(rQG),1)+1);
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȼ<EFBFBD><EFBFBD><EFBFBD>
stateVarCount=2*length(Volt);
measurements=length(mVolt)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG);
fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> %f\n',measurements/stateVarCount);
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>״̬<EFBFBD><EFBFBD><EFBFBD>Ƽ<EFBFBD><EFBFBD><EFBFBD>
% SEVolt=sdpvar(length(Volt),1,'full');
% SEVAngel=sdpvar(length(Vangle),1,'full');
% Objective=(SEVolt-mVolt)'*(1./sigma^2*eye(length(mVolt)))*(SEVolt-mVolt);%<EFBFBD><EFBFBD>ѹ
%% ֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% cmpSEV=SEVolt.*exp(1j*SEVAngel); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
% cmpSEBranchI=(cmpSEV(lineI)-cmpSEV(lineJ))./(lineR+1j*lineX);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% SEBranchI=abs(cmpSEBranchI);% ֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
% Objective=Objective+(SEBranchI-mBranchI)'*(1./sigma^2*eye(length(mBranchI)))*(SEBranchI-mBranchI);%%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% ֧·<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% 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ע<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD>
zerosInjectionIndex=1:length(Volt);
zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
% PQ=diag(cmpSEV)*cmpY*cmpSEV;
% zeroInjP=real(PQ(zerosInjectionIndex));%% 0ע<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD>
% zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0ע<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD>
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ע<EFBFBD>빦<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD>ҵ<EFBFBD>ֻ<EFBFBD>з<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ľڵ<EFBFBD>
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<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% 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 );
opts = optiset('solver','ipopt');
opts.maxiter=3500;
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
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
SEVolt=x(1:length(Volt));
SEVAngel=x(length(Volt)+Balance:2*length(Volt));
fprintf('Ŀ<EFBFBD>꺯<EFBFBD><EFBFBD>Ϊ:%f\n',fval);
fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n');
(abs(rVolt-double(SEVolt)))./(rVolt)
MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)