From cff1dde6f40bd57fee6eeb95e9a77356e9705c56 Mon Sep 17 00:00:00 2001 From: facat Date: Sat, 31 Aug 2013 10:05:49 +0800 Subject: [PATCH] =?UTF-8?q?LM=E6=B3=95=E6=95=88=E6=9E=9C=E4=B9=9F=E4=B8=8D?= =?UTF-8?q?=E5=A5=BD=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- run.m | 114 ++++++++++++++++++++++++++++++++++++++-------------------- 1 file changed, 75 insertions(+), 39 deletions(-) diff --git a/run.m b/run.m index fa012de..823edb0 100644 --- a/run.m +++ b/run.m @@ -72,7 +72,7 @@ onlyQG=setdiff(QGi,PDQDi); %% 计算方差 measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma)); 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)); %% 冗余度计算 stateVarCount=2*length(Volt); @@ -97,7 +97,7 @@ fprintf(' % load('mTransP'); % load('mTransQ'); %% 自己写的微分代码 -% 以下都是Jacobi矩阵 +% 初始化一些值 SEVolt=sparse(ones(length(mVolt),1)); SEVolt(Balance)=rVolt(Balance); SEVAngle=sparse(zeros(length(mVolt),1)); @@ -105,9 +105,13 @@ maxD=1000; Iteration=0; optimalCondition=100; eps=1e-4; +mu=0; +v=2; +ojbFunDecrease=1000;% 目标函数下降 +% 以下都是Jacobi矩阵 % while max(abs(g))>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_dTyta=sparse(length(mVolt),length(mVolt)); @@ -200,33 +204,33 @@ while max(abs(maxD))>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 ; @@ -243,14 +247,19 @@ while max(abs(maxD))>eps z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ]; G=H'*W*H; 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]; + % 利用 Levenber-Marquardt + % if Iteration==0 + % mu=max(diag(G)); + % end a=G; b=-g; - + + % 平衡节点相角恒定 a(length(mVolt)+Balance,:)=0; a(:,length(mVolt)+Balance)=0; @@ -263,14 +272,41 @@ while max(abs(maxD))>eps b(Balance)=0; dX=a\b; 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)))) fprintf('max abs g:%f\n',full(max(abs(g)))); - % 更新变量 SEVolt=SEVolt+dX(1:length(mVolt))*dXStep; - Iteration=Iteration+1; SEVAngle=SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep; 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(Balance)=0; optimalCondition(Balance+length(mVolt))=0; @@ -280,7 +316,7 @@ fprintf(' fval=full((z-h)'*W*(z-h)); fprintf('目标函数为:%f\n',fval); fprintf('相对误差\n'); -(abs(rVolt-double(SEVolt)))./(rVolt) +(abs(rVolt-double(SEVolt)))./(rVolt); MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngle) plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle ) % 检查最优性条件