stateestimation-self-derivatives-measurementdistribution/run.asv

addpath('.\Powerflow')
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0');
%% <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <20><>ѹ <20>ڵ<EFBFBD><DAB5><EFBFBD><EFBFBD><EFBFBD> ֧·<D6A7><C2B7><EFBFBD><EFBFBD> <20>ڵ㹦<DAB5><E3B9A6> ֧·<D6A7><C2B7><EFBFBD><EFBFBD>
%%
%% ״̬<D7B4><CCAC>
% <20><>ѹ <20><><EFBFBD><EFBFBD>
%%
%% <20><>ʼ<EFBFBD><CABC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
sigma=0.05;% <20><>׼<EFBFBD><D7BC>
%% <20><>ѹ
%<25><>ѹ<EFBFBD><D1B9>ֵ
rVolt=Volt; %<25><>ֵ
mVolt=rVolt.*(normrnd(0,sigma,length(Volt),1)+1);%<25><>ѹ<EFBFBD><D1B9><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% <20><><EFBFBD><EFBFBD>
%ע<><D7A2><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
cmpY=Y.*exp(1j*sparse(r,c,Yangle,length(Y),length(Y)));%<25><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɾ<EFBFBD><C9BE><EFBFBD>
cmpV=Volt.*exp(1j*Vangle); %<25><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
cmpI=cmpY*cmpV;% ע<><D7A2><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
rI=abs(cmpI); %ע<><D7A2><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҫ<EFBFBD><D2AA><EFBFBD>ǵ<EFBFBD><C7B5><EFBFBD><EFBFBD><EFBFBD>ֵ
mI=rI.*(normrnd(0,sigma,length(rI),1)+1);%<25><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% ֧·<D6A7><C2B7><EFBFBD><EFBFBD>
% ֧·<D6A7><C2B7><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);%<25><><EFBFBD><EFBFBD>֧·<D6A7><C2B7><EFBFBD><EFBFBD>
rBranchI=abs(cmpBranchI);% ֧·<D6A7><C2B7><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%֧·<D6A7><C2B7><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% ֧·<D6A7><C2B7><EFBFBD><EFBFBD>
rBranchP=real((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI));
mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%֧·<D6A7><C2B7><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);%֧·<D6A7><C2B7><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% ע<>빦<EFBFBD><EBB9A6>
rPD=PD(PD~=0);
PDi=find(PD~=0);
rQD=QD(QD~=0);
QDi=find(QD~=0);
rPG=PG(PG~=0);
PGi=find(PG~=0);
rQG=QG(QG~=0);
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);
%% <20><><EFBFBD><EFBFBD><EFBFBD>ȼ<EFBFBD><C8BC><EFBFBD>
stateVarCount=2*length(Volt);
measurements=length(mVolt)+length(mI)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG);
fprintf('<27><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> %f\n',measurements/stateVarCount);
%% <20><><EFBFBD><EFBFBD>״̬<D7B4><CCAC><EFBFBD>Ƽ<EFBFBD><C6BC><EFBFBD>
SEVolt=sdpvar(length(Volt),1);
SEVAngel=sdpvar(length(Vangle),1);
Objective=(SEVolt-mVolt)'*(1./sigma^2*eye(length(mVolt)))*(SEVolt-mVolt);%%<25><>ѹ
%% ֧·<D6A7><C2B7><EFBFBD><EFBFBD>
cmpSEV=SEVolt.*exp(1j*SEVAngel); %<25><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
cmpSEBranchI=(cmpV(lineI)-cmpV(lineJ))./(lineR+1j*lineX);%<25><><EFBFBD><EFBFBD>֧·<D6A7><C2B7><EFBFBD><EFBFBD>
SEBranchI=abs(cmpSEBranchI);% ֧·<D6A7><C2B7><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
Objective=Objective+(SEBranchI-mBranchI)'*(1./sigma^2*eye(length(mBranchI)))*(SEBranchI-mBranchI);%%<25><><EFBFBD><EFBFBD>
%% ֧·<D6A7><C2B7><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ע<30><D7A2><EFBFBD>ڵ<EFBFBD>
zerosInjectionIndex=1:length(Volt);
zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
PQ=diag(SEVolt)*conj(cmpY*SEVolt);
zeroInjP=real(PQ(zerosInjectionIndex));%% 0ע<30><D7A2><EFBFBD>ڵ<EFBFBD>
zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0ע<30><D7A2><EFBFBD>ڵ<EFBFBD>
%% <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ע<EFBFBD>빦<EFBFBD><EBB9A6>
% <20><><EFBFBD>ҵ<EFBFBD>ֻ<EFBFBD>з<EFBFBD><D0B7><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ľڵ<C4BD>
PDQDi=union(PDi,QDi);
onlyPG=setdiff(PGi,PDQDi);
onlyQG=setdiff(QGi,PDQDi);
PG-real(PQ())
%% YALMIP<49><50><EFBFBD><EFBFBD>
Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0];
% Constraints=[[zeros(length(c)) A' -eye(length(lbounds))]*x==-c;[A zeros(length(b)) zeros(length(b),length(lbounds))]*x<=b;10>=x>=0];
options = sdpsettings('verbose',2,'solver','ipopt');
solvesdp(Constraints,Objective,options)
double(Objective)
fprintf('<27><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n');
(abs(rVolt-double(SEVolt)))./(rVolt)