修改了权重估计结果很好，但是对于1047节点，收敛性不好。准备用LM法。
Signed-off-by: facat <facat@facat.cn>
This commit is contained in:
facat
parent
4ba940dc0a
commit
62389acec5
|
|
@ -4,9 +4,11 @@ subplot(2,1,1);
|
|||
plot(x,rVolt,'r-');
|
||||
hold on
|
||||
plot(x,SEVolt);
|
||||
legend('real','se');
|
||||
subplot(2,1,2);
|
||||
plot(x,rVAngel,'r-');
|
||||
hold on
|
||||
plot(x,SEVAngel);
|
||||
legend('real','se');
|
||||
end
|
||||
|
||||
|
|
|
|||
76
run.m
76
run.m
|
|
@ -4,7 +4,7 @@ clc
|
|||
addpath('.\Powerflow')
|
||||
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('s1047.dat', '0');
|
||||
%% 开始生成量测量
|
||||
sigma=0.05;% 标准差
|
||||
sigma=0.03;% 标准差
|
||||
%% 电压
|
||||
%电压幅值
|
||||
rVolt=Volt; %幅值
|
||||
|
|
@ -70,9 +70,9 @@ PDQDi=union(PDi,QDi);
|
|||
onlyPG=setdiff(PGi,PDQDi);
|
||||
onlyQG=setdiff(QGi,PDQDi);
|
||||
%% 计算方差
|
||||
measureSigma=([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma).^2;
|
||||
measureSigma(measureSigma<1e-4)=mean(measureSigma);
|
||||
W=diag(1./measureSigma) ;
|
||||
measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma));
|
||||
measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6));
|
||||
W=diag(1./measureSigma.^2) ;
|
||||
% W=sparse(1:length(W),1:length(W),400,length(W),length(W));
|
||||
%% 冗余度计算
|
||||
stateVarCount=2*length(Volt);
|
||||
|
|
@ -107,7 +107,7 @@ optimalCondition=100;
|
|||
eps=1e-4;
|
||||
% 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));
|
||||
|
|
@ -200,33 +200,33 @@ while max(abs(optimalCondition))>eps
|
|||
,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);
|
||||
% 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);
|
||||
%% 进入迭代
|
||||
H=[dV_dV,dV_dTyta;
|
||||
dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj ;
|
||||
|
|
@ -245,9 +245,12 @@ while max(abs(optimalCondition))>eps
|
|||
g=-H'*W*(z-h);
|
||||
|
||||
% 形成大的求解矩阵
|
||||
a=[G C';
|
||||
C zeros(size(C,1),size(C,1))];
|
||||
b=[-g;-c];
|
||||
% a=[G C';
|
||||
% C zeros(size(C,1),size(C,1))];
|
||||
% b=[-g;-c];
|
||||
a=G;
|
||||
b=-g;
|
||||
|
||||
% 平衡节点相角恒定
|
||||
a(length(mVolt)+Balance,:)=0;
|
||||
a(:,length(mVolt)+Balance)=0;
|
||||
|
|
@ -268,8 +271,7 @@ while max(abs(optimalCondition))>eps
|
|||
Iteration=Iteration+1;
|
||||
SEVAngle=SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep;
|
||||
lamda=-dX(length(mVolt)*2+1:end);
|
||||
optimalCondition=[-g+C'*lamda;
|
||||
c];
|
||||
optimalCondition=-g;
|
||||
optimalCondition(Balance)=0;
|
||||
optimalCondition(Balance+length(mVolt))=0;
|
||||
end
|
||||
|
|
|
|||
Loading…
Reference in New Issue