迭代很多次才能满足最优性条件，干脆直接判断最大修正量。

Signed-off-by: dmy@lab <dmy@lab.lab>

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 ;
+        -dPdV(PDi,:),-dPdTyta(PDi,:);
-        dLQij_dVi+dLQij_dVj,dLQij_dThetai+dLQij_dThetaj ;
+        -dQdV(QDi,:),-dQdTyta(QDi,:)];%jacobi
-        dTPij_dVi+dTPij_dVj,dTPij_dThetai+dTPij_dThetaj;
-        dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj];%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];
+    SEPD=-diag(SEVolt)*Y.* ( spfun (@cos, VAngleIJ ) )*SEVolt;
-    z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ];
+    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;

公式/公式.tex

@@ -201,7 +201,22 @@ Q_{ij}&=-\frac{V_1^2}{k^2}B_{ij}-\frac{V_1}{k} V_2[sin(\theta_1 - \theta_2)G_{ij
 \end{aligned}
 \end{equation}
 
-ps.已检验过线路的公式。
+为了推潮流公式，先从简单的开始推起。
+\begin{array}{c}
+	1 \\
+	d \\
+\end{array}
+
+潮流方程有功的公式为
+\begin{equation}
+\Delta P=diag(V)[Y.*cos(\theta e^T -e \theta^T-\alpha)]V+P_{D}-P_{G}=0
+\end{equation}
+
+\begin{equation}
+\frac{\partial P}{\partial V}=
+\end{equation}
+
+ps.已检验过线路以及变压器支路的公式。
 \par
 以上公式已经可以完成状态估计，若要实现更好的收敛性，需要利用二阶导数。
 \par