加入Armijo搜索,效果不好。

Signed-off-by: facat <facat@facat.cn>

Armijo.m

@@ -0,0 +1,17 @@
+function [ output_args ] = Armijo(z,newwordParameter,W,SEVolt,SEVAngle,d,g )
+%ArmijoËÑË÷
+
+a=0.9;
+b=0.1;
+output_args=a;
+for m=1:100
+    objfun( SEVolt,SEVAngle,W,z,newwordParameter )
+    if objfun( SEVolt+a^m*d(1:length(SEVolt)),SEVAngle+a^m*d(length(SEVolt)+1:end),W,z,newwordParameter )<=objfun( SEVolt,SEVAngle,W,z,newwordParameter )+b*a^m*g'*d
+        output_args= a^m;
+        return;
+    end
+end
+
+
+end
+

objfun.m

@@ -0,0 +1,18 @@
+function [ output_args ] = objfun( SEVolt,SEVAngle,W,z,newwordParameter )
+%Ä¿±êº¯Êý
+%%
+lineI=newwordParameter.line.lineI;
+lineJ=newwordParameter.line.lineJ;
+lineR=newwordParameter.line.lineR;
+lineX=newwordParameter.line.lineX;
+lineB2=newwordParameter.line.lineB2;
+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];
+output_args=(z-h)'*W*(z-h);
+end
+

run.m

@@ -2,15 +2,9 @@ clear
 clc
 % yalmip('clear')
 addpath('.\Powerflow')
-[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee4.dat', '0');
+[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('S1047.dat', '0');
-%% 量测量
-% 电压 节点电流 支路电流 节点功率 支路功率
-%%
-%% 状态量
-% 电压 相角
-%%
 %% 开始生成量测量
-sigma=0.05;% 标准差
+sigma=0.09;% 标准差
 %% 电压
 %电压幅值
 rVolt=Volt; %幅值
@@ -41,6 +35,18 @@ rBranchP= BranchP( cmpV,cmpBranchI,lineI,lineB2 );
 mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%支路功率量测量
 rBranchQ=BranchQ( cmpV,cmpBranchI,lineI,lineB2 );
 mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量
+%% 变压器功率
+transI=newwordParameter.trans.transI;
+transJ=newwordParameter.trans.transJ;
+transK=newwordParameter.trans.transK;
+transR=newwordParameter.trans.transR;
+transX=newwordParameter.trans.transX;
+transG=real(1./(transR+1j*transX));
+transB=imag(1./(transR+1j*transX));
+rTransP=TransPower( newwordParameter,rVolt,rVAngel );
+rTransQ=TransReactivePower( newwordParameter,rVolt,rVAngel );
+mTransP=rTransP.*(normrnd(0,sigma,length(rTransP),1)+1);
+mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1);
 %% 注入功率
 rPD=PD;
 PDi=find(PD~=0);
@@ -57,29 +63,18 @@ mQG=rQG.*(normrnd(0,sigma,length(rQG),1)+1);
 %% 0注入节点
 zerosInjectionIndex=1:length(Volt);
 zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
-%% 变压器功率
-transI=newwordParameter.trans.transI;
-transJ=newwordParameter.trans.transJ;
-transK=newwordParameter.trans.transK;
-transR=newwordParameter.trans.transR;
-transX=newwordParameter.trans.transX;
-transG=real(1./(transR+1j*transX));
-transB=imag(1./(transR+1j*transX));
-rTransP=TransPower( newwordParameter,rVolt,rVAngel );
-rTransQ=TransReactivePower( newwordParameter,rVolt,rVAngel );
-mTransP=rTransP.*(normrnd(0,sigma,length(rTransP),1)+1);
-mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1);
 %% 发电机注入功率
 % 先找到只有发电机的节点
 PDQDi=union(PDi,QDi);
 onlyPG=setdiff(PGi,PDQDi);
 onlyQG=setdiff(QGi,PDQDi);
+%% 计算方差
+W=(diag([rVolt;rBranchP;rBranchQ;rTransP;rTransQ])*sigma).^2;
+W=sparse(1:length(W),1:length(W),400,length(W),length(W));
 %% 冗余度计算
 stateVarCount=2*length(Volt);
 measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG)+length(mTransP)+length(mTransQ);
 fprintf('冗余度 %f\n',measurements/stateVarCount);
-%% Opti ToolBox
-Busnum=length(Volt);
 %% save
 % save('mVolt','mVolt');
 % save('mPG','mPG');
@@ -101,118 +96,116 @@ Busnum=length(Volt);
 %% 自己写的微分代码
 % 以下都是Jacobi矩阵
 SEVolt=sparse(ones(length(mVolt),1));
-SEVAngel=sparse(zeros(length(mVolt),1));
+SEVAngle=sparse(zeros(length(mVolt),1));
 maxD=1000;
 Iteration=0;
-while maxD>1e-3;
+g=100;
+while max(abs(g))>1e-5;
+% while g>1e-5
     % 电压
     dV_dV=sparse(1:length(mVolt),1:length(mVolt),1,length(mVolt),length(mVolt));%电压量测量的微分
     dV_dTyta=sparse(length(mVolt),length(mVolt));
     % 线路支路
     dLPij_dVi=sparse(1:length(lineI),lineI, ...
         -SEVolt(lineJ).*( ...
-        lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
+        lineG.*cos(SEVAngle(lineI)-SEVAngle(lineJ)) +lineB.*sin(SEVAngle(lineI)-SEVAngle(lineJ))...
         )...
         +2*(lineG).*SEVolt(lineI) ...
         ,length(lineI),length(mVolt));%线路的
     dLPij_dVj=sparse(1:length(lineI),lineJ, ...
         -SEVolt(lineI).*( ...
-        lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
+        lineG.*cos(SEVAngle(lineI)-SEVAngle(lineJ)) +lineB.*sin(SEVAngle(lineI)-SEVAngle(lineJ))...
         ) ...
         ,length(lineI),length(mVolt));%线路的
     dLPij_dThetai=sparse(1:length(lineI),lineI, ...
         SEVolt(lineI).*SEVolt(lineJ).*( ...
-        lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
+        lineG.*sin(SEVAngle(lineI)-SEVAngle(lineJ)) -lineB.*cos(SEVAngle(lineI)-SEVAngle(lineJ))...
         )...
         ,length(lineI),length(mVolt));%线路的
     dLPij_dThetaj=sparse(1:length(lineI),lineJ, ...
         -SEVolt(lineI).*SEVolt(lineJ).*( ...
-        lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
+        lineG.*sin(SEVAngle(lineI)-SEVAngle(lineJ)) -lineB.*cos(SEVAngle(lineI)-SEVAngle(lineJ))...
         )...
         ,length(lineI),length(mVolt));%线路的
     
     dLQij_dVi=sparse(1:length(lineI),lineI, ...
         -SEVolt(lineJ).*( ...
-        lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
+        lineG.*sin(SEVAngle(lineI)-SEVAngle(lineJ)) -lineB.*cos(SEVAngle(lineI)-SEVAngle(lineJ))...
         )...
         -2*(lineB+lineB2).*SEVolt(lineI) ...
         ,length(lineI),length(mVolt));%线路的
     dLQij_dVj=sparse(1:length(lineI),lineJ, ...
         -SEVolt(lineI).*( ...
-        lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
+        lineG.*sin(SEVAngle(lineI)-SEVAngle(lineJ)) -lineB.*cos(SEVAngle(lineI)-SEVAngle(lineJ))...
         ) ...
         ,length(lineI),length(mVolt));%线路的
     dLQij_dThetai=sparse(1:length(lineI),lineI, ...
         -SEVolt(lineI).*SEVolt(lineJ).*( ...
-        lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
+        lineG.*cos(SEVAngle(lineI)-SEVAngle(lineJ)) +lineB.*sin(SEVAngle(lineI)-SEVAngle(lineJ))...
         ) ...
         ,length(lineI),length(mVolt));%线路的
     dLQij_dThetaj=sparse(1:length(lineI),lineJ, ...
         SEVolt(lineI).*SEVolt(lineJ).*( ...
-        lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
+        lineG.*cos(SEVAngle(lineI)-SEVAngle(lineJ)) +lineB.*sin(SEVAngle(lineI)-SEVAngle(lineJ))...
         ) ...
         ,length(lineI),length(mVolt));%线路的
     % 变压器支路
     dTPij_dVi=sparse(1:length(transI),transI, ...
         -SEVolt(transJ)./transK.*( ...
-        transG.*cos(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))...
+        transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))...
         )...
         +2*transG.*SEVolt(transI)./(transK.^2) ...
         ,length(transI),length(mVolt));%变压器
     dTPij_dVj=sparse(1:length(transI),transJ, ...
         -SEVolt(transI)./transK.*( ...
-        transG.*cos(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))...
+        transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))...
         ) ...
         ,length(transI),length(mVolt));%变压器
     dTPij_dThetai=sparse(1:length(transI),transI, ...
         SEVolt(transI)./transK.*SEVolt(transJ).*( ...
-        transG.*sin(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(transJ))...
+        transG.*sin(SEVAngle(transI)-SEVAngle(transJ)) -transB.*cos(SEVAngle(transI)-SEVAngle(transJ))...
         )...
         ,length(transI),length(mVolt));%变压器
     dTPij_dThetaj=sparse(1:length(transI),transJ, ...
         -SEVolt(transI)./transK.*SEVolt(transJ).*( ...
-        transG.*sin(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(transJ))...
+        transG.*sin(SEVAngle(transI)-SEVAngle(transJ)) -transB.*cos(SEVAngle(transI)-SEVAngle(transJ))...
         )...
         ,length(transI),length(mVolt));%变压器
     dTQij_dVi=sparse(1:length(transI),transI, ...
         -SEVolt(transJ)./transK.*( ...
-        transG.*sin(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(transJ))...
+        transG.*sin(SEVAngle(transI)-SEVAngle(transJ)) -transB.*cos(SEVAngle(transI)-SEVAngle(transJ))...
         )...
         -2*transB.*SEVolt(transI)./(transK.^2) ...
         ,length(transI),length(mVolt));%变压器
     dTQij_dVj=sparse(1:length(transI),transJ, ...
         -SEVolt(transI)./transK.*( ...
-        transG.*sin(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(transJ))...
+        transG.*sin(SEVAngle(transI)-SEVAngle(transJ)) -transB.*cos(SEVAngle(transI)-SEVAngle(transJ))...
         ) ...
         ,length(transI),length(mVolt));%变压器
     dTQij_dThetai=sparse(1:length(transI),transI, ...
         -SEVolt(transI)./transK.*SEVolt(transJ).*( ...
-        transG.*cos(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))...
+        transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))...
         ) ...
         ,length(transI),length(mVolt));%变压器
     dTQij_dThetaj=sparse(1:length(transI),transJ, ...
         SEVolt(transI)./transK.*SEVolt(transJ).*( ...
-        transG.*cos(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))...
+        transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))...
         ) ...
         ,length(transI),length(mVolt));%变压器
     %% 进入迭代
-    % H=[dV_dV ,dV_dTyta;
-    %   dLPij_dVi+dLPij_dVj, dLPij_dThetai+dLPij_dThetaj;
-    %   dLQij_dVi+dLQij_dVj, dLQij_dThetai+dLQij_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
     
-    SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngel),lineI,lineJ,lineR,lineX );%复数支路电流
+    SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngle),lineI,lineJ,lineR,lineX );%复数支路电流
-    SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngel),SEBranchI,lineI,lineB2 );
+    SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
-    SEBranchQ=BranchQ( SEVolt.*exp(1j*SEVAngel),SEBranchI,lineI,lineB2 );
+    SEBranchQ=BranchQ( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
-    SETransP=TransPower( newwordParameter,SEVolt,SEVAngel );
+    SETransP=TransPower( newwordParameter,SEVolt,SEVAngle );
-    SETransQ=TransReactivePower( newwordParameter,SEVolt,SEVAngel );
+    SETransQ=TransReactivePower( newwordParameter,SEVolt,SEVAngle );
     h=[SEVolt;SEBranchP;SEBranchQ;SETransP;SETransQ];
     z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ];
-    W=sparse(1:length(h),1:length(h),1/sigma.^2,length(h),length(h));
     G=H'*W*H;
     g=-H'*W*(z-h);
     % 平衡节点相角恒定;
@@ -221,11 +214,13 @@ while maxD>1e-3;
     G=G+sparse(length(mVolt)+Balance,length(mVolt)+Balance,1,length(mVolt)*2,length(mVolt)*2);
     g(length(mVolt)+Balance)=0;
     dX=G\-g;
+    dXStep=1;
+%     dXStep=Armijo(z,newwordParameter,W,SEVolt,SEVAngle,dX,g );
     maxD=max(abs(dX))
     % 更新变量
-    SEVolt=SEVolt+dX(1:length(mVolt));
+    SEVolt=SEVolt+dX(1:length(mVolt))*dXStep;
     Iteration=Iteration+1;
-    SEVAngel=SEVAngel+dX(length(mVolt)+1:end);
+    SEVAngle=SEVAngle+dX(length(mVolt)+1:end)*dXStep;
 end
 %% 输出结果
 fprintf('迭代%d次\n',Iteration);
@@ -233,5 +228,11 @@ fval=full((z-h)'*W*(z-h));
 fprintf('目标函数为:%f\n',fval);
 fprintf('相对误差\n');
 (abs(rVolt-double(SEVolt)))./(rVolt)
-MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
+MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngle)
-plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngel )
+plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle )
+% 检查最优性条件
+if any(abs(g)>1e-5)
+    fprintf('最优性条件不满足\n')
+else
+    fprintf('最优性条件满足\n')
+end