Merge branch 'master' of http://emp2:3000/facat/stateestimation-self-derivatives.git
Conflicts: 鍏?寮?/鍏?寮?.tex
This commit is contained in:
dmy@lab
commit
df1943a96d
65
run.m
65
run.m
|
|
@ -2,7 +2,7 @@ clear
|
|||
clc
|
||||
% yalmip('clear')
|
||||
addpath('.\Powerflow')
|
||||
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee118.dat', '0');
|
||||
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('E:\算例\feeder33\feeder33ieee.txt', '0');
|
||||
%% 开始生成量测量
|
||||
sigma=0.03;% 标准差
|
||||
%% 电压
|
||||
|
|
@ -70,7 +70,8 @@ PDQDi=union(PDi,QDi);
|
|||
onlyPG=setdiff(PGi,PDQDi);
|
||||
onlyQG=setdiff(QGi,PDQDi);
|
||||
%% 计算方差
|
||||
measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma));
|
||||
% measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma));
|
||||
measureSigma=abs(([rVolt;rPD(PDi);rQD(QDi);].*sigma));
|
||||
measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6));
|
||||
W=sparse(diag(1./measureSigma.^2)) ;
|
||||
% W=sparse(1:length(W),1:length(W),400,length(W),length(W));
|
||||
|
|
@ -100,7 +101,7 @@ fprintf('
|
|||
% 初始化一些值
|
||||
SEVolt=sparse(ones(length(mVolt),1));
|
||||
SEVolt(Balance)=rVolt(Balance);
|
||||
SEVAngle=sparse(-0.02*ones(length(mVolt),1));
|
||||
SEVAngle=sparse(-0.00*ones(length(mVolt),1));
|
||||
maxD=1000;
|
||||
Iteration=0;
|
||||
optimalCondition=100;
|
||||
|
|
@ -111,7 +112,7 @@ ojbFunDecrease=1000;% Ŀ
|
|||
% 以下都是Jacobi矩阵
|
||||
% while max(abs(g))>1e-5;
|
||||
% while maxD>1e-5
|
||||
while max(abs(optimalCondition))>eps
|
||||
while max(abs(maxD))>eps
|
||||
% 电压
|
||||
dV_dV=sparse(1:length(mVolt),1:length(mVolt),1,length(mVolt),length(mVolt));%电压量测量的微分
|
||||
dV_dTyta=sparse(length(mVolt),length(mVolt));
|
||||
|
|
@ -202,6 +203,36 @@ while max(abs(optimalCondition))>eps
|
|||
transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))...
|
||||
) ...
|
||||
,length(transI),length(mVolt));%变压器
|
||||
%% 考虑注入功率
|
||||
% 等式约束的Jacobi
|
||||
r=newwordParameter.r;
|
||||
c=newwordParameter.c;
|
||||
Yangle=newwordParameter.Yangle;
|
||||
VAngleIJ=sparse(r,c,SEVAngle(r)-SEVAngle(c) -Yangle,length(mVolt),length(mVolt)) ;
|
||||
YdotSin=Y.* ( spfun(@sin,VAngleIJ) );
|
||||
YdotCos=Y.* ( spfun (@cos, VAngleIJ ) );
|
||||
diag_Volt_YdotCos=diag(SEVolt)*YdotCos;
|
||||
diag_Volt_YdotSin=diag(SEVolt)*YdotSin;
|
||||
YdotCosVolt=YdotCos*Volt;
|
||||
YdotSinVolt=YdotSin*Volt;
|
||||
diag_Volt_YdotCosVolt=diag_Volt_YdotCos*Volt;
|
||||
diag_Volt_YdotSinVolt=diag_Volt_YdotSin*Volt;
|
||||
diag_YdotSinVolt_=diag(YdotSinVolt);
|
||||
diag_YdotCosVolt_=diag(YdotCosVolt);
|
||||
dPdTyta=diag_Volt_YdotSin*diag(SEVolt)-diag_YdotSinVolt_*diag(SEVolt); % 简化第三次
|
||||
dQdTyta=-diag_Volt_YdotCos*diag(SEVolt)+diag_YdotCosVolt_*diag(SEVolt);%dQ/dThyta
|
||||
dPdV=diag_YdotCosVolt_+diag_Volt_YdotCos;%dP/dV
|
||||
dQdV=diag_YdotSinVolt_+diag_Volt_YdotSin;%dQ/dV
|
||||
% % C 是等式约束 c 的Jacobi
|
||||
% C=[dPdV dPdTyta;
|
||||
% dQdV dQdTyta];
|
||||
% C=C(zerosInjectionIndex,:);
|
||||
% % 形成等式约束 c
|
||||
% nodeP=diag_Volt_YdotCosVolt;
|
||||
% nodeQ=diag_Volt_YdotSinVolt;
|
||||
% nodePQ=[nodeP;nodeQ];
|
||||
% c=nodePQ(zerosInjectionIndex);
|
||||
|
||||
%% 考虑等式约束
|
||||
% 等式约束的Jacobi
|
||||
% r=newwordParameter.r;
|
||||
|
|
@ -232,19 +263,31 @@ while max(abs(optimalCondition))>eps
|
|||
% nodePQ=[nodeP;nodeQ];
|
||||
% c=nodePQ(zerosInjectionIndex);
|
||||
%% 进入迭代
|
||||
% H=[dV_dV,dV_dTyta;
|
||||
% dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj ;
|
||||
% dLQij_dVi+dLQij_dVj,dLQij_dThetai+dLQij_dThetaj ;
|
||||
% dTPij_dVi+dTPij_dVj,dTPij_dThetai+dTPij_dThetaj;
|
||||
% dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj];%jacobi
|
||||
H=[dV_dV,dV_dTyta;
|
||||
dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj ;
|
||||
dLQij_dVi+dLQij_dVj,dLQij_dThetai+dLQij_dThetaj ;
|
||||
dTPij_dVi+dTPij_dVj,dTPij_dThetai+dTPij_dThetaj;
|
||||
dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj];%jacobi
|
||||
-dPdV(PDi,:),-dPdTyta(PDi,:);
|
||||
-dQdV(QDi,:),-dQdTyta(QDi,:)];%jacobi
|
||||
|
||||
SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngle),lineI,lineJ,lineR,lineX );%复数支路电流
|
||||
SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
|
||||
SEBranchQ=BranchQ( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
|
||||
SETransP=TransPower( newwordParameter,SEVolt,SEVAngle );
|
||||
SETransQ=TransReactivePower( newwordParameter,SEVolt,SEVAngle );
|
||||
h=[SEVolt;SEBranchP;SEBranchQ;SETransP;SETransQ];
|
||||
z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ];
|
||||
SEPD=-diag(SEVolt)*Y.* ( spfun (@cos, VAngleIJ ) )*SEVolt;
|
||||
SEQD=-diag(SEVolt)*Y.* ( spfun(@sin,VAngleIJ) )*SEVolt;
|
||||
|
||||
h=[SEVolt;SEPD(PDi);SEQD(QDi);];
|
||||
|
||||
% h=[SEVolt;SEBranchP;SEBranchQ;SETransP;SETransQ];
|
||||
% z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ];
|
||||
|
||||
z=[mVolt;mPD(PDi);mQD(QDi)];
|
||||
|
||||
|
||||
G=H'*W*H;
|
||||
g=-H'*W*(z-h);
|
||||
|
||||
|
|
@ -273,7 +316,7 @@ while max(abs(optimalCondition))>eps
|
|||
dX=a\b;
|
||||
dXStep=1;
|
||||
% dXStep=Armijo(z,newwordParameter,W,SEVolt,SEVAngle,dX,g );
|
||||
maxD=max(abs(dX(1:length(mVolt))))
|
||||
maxD=max(abs(dX))
|
||||
fprintf('max abs g:%f\n',full(max(abs(g))));
|
||||
SEVolt=SEVolt+dX(1:length(mVolt))*dXStep;
|
||||
SEVAngle=SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep;
|
||||
|
|
|
|||
Loading…
Reference in New Issue