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

LM法效果也不好。 

Signed-off-by: facat <facat@facat.cn>
This commit is contained in:
facat 2013-08-31 10:05:49 +08:00
parent 62389acec5
commit cff1dde6f4
1 changed files with 75 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

110
run.m
View File

@ -72,7 +72,7 @@ onlyQG=setdiff(QGi,PDQDi);
%% 计算方差 %% 计算方差
measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma)); measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma));
measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6)); measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6));
W=diag(1./measureSigma.^2) ; W=sparse(diag(1./measureSigma.^2)) ;
% W=sparse(1:length(W),1:length(W),400,length(W),length(W)); % W=sparse(1:length(W),1:length(W),400,length(W),length(W));
%% 冗余度计算 %% 冗余度计算
stateVarCount=2*length(Volt); stateVarCount=2*length(Volt);
@ -97,7 +97,7 @@ fprintf('
% load('mTransP'); % load('mTransP');
% load('mTransQ'); % load('mTransQ');
%% 自己写的微分代码 %% 自己写的微分代码
% 以下都是Jacobi矩阵 % 初始化一些值
SEVolt=sparse(ones(length(mVolt),1)); SEVolt=sparse(ones(length(mVolt),1));
SEVolt(Balance)=rVolt(Balance); SEVolt(Balance)=rVolt(Balance);
SEVAngle=sparse(zeros(length(mVolt),1)); SEVAngle=sparse(zeros(length(mVolt),1));
@ -105,9 +105,13 @@ maxD=1000;
Iteration=0; Iteration=0;
optimalCondition=100; optimalCondition=100;
eps=1e-4; eps=1e-4;
mu=0;
v=2;
ojbFunDecrease=1000;% 目标函数下降
% 以下都是Jacobi矩阵
% while max(abs(g))>1e-5; % while max(abs(g))>1e-5;
% while maxD>1e-5 % while maxD>1e-5
while max(abs(maxD))>eps while max(abs(optimalCondition))>eps
% 电压 % 电压
dV_dV=sparse(1:length(mVolt),1:length(mVolt),1,length(mVolt),length(mVolt));%电压量测量的微分 dV_dV=sparse(1:length(mVolt),1:length(mVolt),1,length(mVolt),length(mVolt));%电压量测量的微分
dV_dTyta=sparse(length(mVolt),length(mVolt)); dV_dTyta=sparse(length(mVolt),length(mVolt));
@ -200,33 +204,33 @@ while max(abs(maxD))>eps
,length(transI),length(mVolt));%变压器 ,length(transI),length(mVolt));%变压器
%% 考虑等式约束 %% 考虑等式约束
% 等式约束的Jacobi % 等式约束的Jacobi
% r=newwordParameter.r; % r=newwordParameter.r;
% c=newwordParameter.c; % c=newwordParameter.c;
% Yangle=newwordParameter.Yangle; % Yangle=newwordParameter.Yangle;
% VAngleIJ=sparse(r,c,SEVAngle(r)-SEVAngle(c) -Yangle,length(mVolt),length(mVolt)) ; % VAngleIJ=sparse(r,c,SEVAngle(r)-SEVAngle(c) -Yangle,length(mVolt),length(mVolt)) ;
% YdotSin=Y.* ( spfun(@sin,VAngleIJ) ); % YdotSin=Y.* ( spfun(@sin,VAngleIJ) );
% YdotCos=Y.* ( spfun (@cos, VAngleIJ ) ); % YdotCos=Y.* ( spfun (@cos, VAngleIJ ) );
% diag_Volt_YdotCos=diag(SEVolt)*YdotCos; % diag_Volt_YdotCos=diag(SEVolt)*YdotCos;
% diag_Volt_YdotSin=diag(SEVolt)*YdotSin; % diag_Volt_YdotSin=diag(SEVolt)*YdotSin;
% YdotCosVolt=YdotCos*Volt; % YdotCosVolt=YdotCos*Volt;
% YdotSinVolt=YdotSin*Volt; % YdotSinVolt=YdotSin*Volt;
% diag_Volt_YdotCosVolt=diag_Volt_YdotCos*Volt; % diag_Volt_YdotCosVolt=diag_Volt_YdotCos*Volt;
% diag_Volt_YdotSinVolt=diag_Volt_YdotSin*Volt; % diag_Volt_YdotSinVolt=diag_Volt_YdotSin*Volt;
% diag_YdotSinVolt_=diag(YdotSinVolt); % diag_YdotSinVolt_=diag(YdotSinVolt);
% diag_YdotCosVolt_=diag(YdotCosVolt); % diag_YdotCosVolt_=diag(YdotCosVolt);
% dPdTyta=diag_Volt_YdotSin*diag(SEVolt)-diag_YdotSinVolt_*diag(SEVolt); % 简化第三次 % dPdTyta=diag_Volt_YdotSin*diag(SEVolt)-diag_YdotSinVolt_*diag(SEVolt); % 简化第三次
% dQdTyta=-diag_Volt_YdotCos*diag(SEVolt)+diag_YdotCosVolt_*diag(SEVolt);%dQ/dThyta % dQdTyta=-diag_Volt_YdotCos*diag(SEVolt)+diag_YdotCosVolt_*diag(SEVolt);%dQ/dThyta
% dPdV=diag_YdotCosVolt_+diag_Volt_YdotCos;%dP/dV % dPdV=diag_YdotCosVolt_+diag_Volt_YdotCos;%dP/dV
% dQdV=diag_YdotSinVolt_+diag_Volt_YdotSin;%dQ/dV % dQdV=diag_YdotSinVolt_+diag_Volt_YdotSin;%dQ/dV
% % C 是等式约束 c 的Jacobi % % C 是等式约束 c 的Jacobi
% C=[dPdV dPdTyta; % C=[dPdV dPdTyta;
% dQdV dQdTyta]; % dQdV dQdTyta];
% C=C(zerosInjectionIndex,:); % C=C(zerosInjectionIndex,:);
% % 形成等式约束 c % % 形成等式约束 c
% nodeP=diag_Volt_YdotCosVolt; % nodeP=diag_Volt_YdotCosVolt;
% nodeQ=diag_Volt_YdotSinVolt; % nodeQ=diag_Volt_YdotSinVolt;
% nodePQ=[nodeP;nodeQ]; % nodePQ=[nodeP;nodeQ];
% c=nodePQ(zerosInjectionIndex); % c=nodePQ(zerosInjectionIndex);
%% 进入迭代 %% 进入迭代
H=[dV_dV,dV_dTyta; H=[dV_dV,dV_dTyta;
dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj ; dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj ;
@ -245,12 +249,17 @@ while max(abs(maxD))>eps
g=-H'*W*(z-h); g=-H'*W*(z-h);
% 形成大的求解矩阵 % 形成大的求解矩阵
% a=[G C'; % a=[G C';
% C zeros(size(C,1),size(C,1))]; % C zeros(size(C,1),size(C,1))];
% b=[-g;-c]; % b=[-g;-c];
% 利用 Levenber-Marquardt
% if Iteration==0
% mu=max(diag(G));
% end
a=G; a=G;
b=-g; b=-g;
% 平衡节点相角恒定 % 平衡节点相角恒定
a(length(mVolt)+Balance,:)=0; a(length(mVolt)+Balance,:)=0;
a(:,length(mVolt)+Balance)=0; a(:,length(mVolt)+Balance)=0;
@ -263,14 +272,41 @@ while max(abs(maxD))>eps
b(Balance)=0; b(Balance)=0;
dX=a\b; dX=a\b;
dXStep=1; dXStep=1;
% dXStep=Armijo(z,newwordParameter,W,SEVolt,SEVAngle,dX,g ); % dXStep=Armijo(z,newwordParameter,W,SEVolt,SEVAngle,dX,g );
maxD=max(abs(dX(1:length(mVolt)))) maxD=max(abs(dX(1:length(mVolt))))
fprintf('max abs g:%f\n',full(max(abs(g)))); fprintf('max abs g:%f\n',full(max(abs(g))));
% 更新变量
SEVolt=SEVolt+dX(1:length(mVolt))*dXStep; SEVolt=SEVolt+dX(1:length(mVolt))*dXStep;
Iteration=Iteration+1;
SEVAngle=SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep; SEVAngle=SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep;
lamda=-dX(length(mVolt)*2+1:end); lamda=-dX(length(mVolt)*2+1:end);
% % 计算目标函数下降量
% preObjFun=(z-h)'*W*(z-h);
% oldSEVolt=SEVolt;
% oldSEVAngle=SEVAngle;
% 求更新控制指标 p
% p=(objfun( SEVolt,SEVAngle,W,z,newwordParameter ) - objfun( SEVolt+dX(1:length(mVolt))*dXStep,SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep,W,z,newwordParameter ) )/( .5*dX'*(mu*dX-g) );
% if p>0 %接受更新
% % 更新变量
% SEVolt=SEVolt+dX(1:length(mVolt))*dXStep;
% SEVAngle=SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep;
% lamda=-dX(length(mVolt)*2+1:end);
% mu=mu*max([1/3,1-(2*p-1)^3]);
% v=2;
% else
% mu=mu*v;
% v=2*v;
% end
Iteration=Iteration+1;
% % 更新后目标函数
% h=[SEVolt;SEBranchP;SEBranchQ;SETransP;SETransQ];
% ojbFunDecrease=preObjFun-(z-h)'*W*(z-h)
% if ojbFunDecrease<1e-3
% mu=100000;
% else
% mu=0;
% end
optimalCondition=-g; optimalCondition=-g;
optimalCondition(Balance)=0; optimalCondition(Balance)=0;
optimalCondition(Balance+length(mVolt))=0; optimalCondition(Balance+length(mVolt))=0;
@ -280,7 +316,7 @@ fprintf('
fval=full((z-h)'*W*(z-h)); fval=full((z-h)'*W*(z-h));
fprintf('目标函数为:%f\n',fval); fprintf('目标函数为:%f\n',fval);
fprintf('相对误差\n'); fprintf('相对误差\n');
(abs(rVolt-double(SEVolt)))./(rVolt) (abs(rVolt-double(SEVolt)))./(rVolt);
MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngle) MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngle)
plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle ) plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle )
% 检查最优性条件 % 检查最优性条件