更新电流，有功，无功计算公式。开始考虑接地支路。

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

BranchI.m

@@ -1,5 +1,5 @@
-function [ output_args ] = BranchI( V,lineI,lineJ,lineR,lineX )
+function [ output_args ] = BranchI( V,lineI,lineJ,lineR,lineX,lineB2 )
-output_args=(V(lineI)-V(lineJ))./(lineR+1j*lineX);
+output_args=(V(lineI)-V(lineJ))./(lineR+1j*lineX)+V(lineI).*lineB2;
 
 end
 

BranchP.m

@@ -1,4 +1,4 @@
 function [ output_args ] = BranchP( V,I,lineI,lineB2 )
-output_args=real((V(lineI)).*conj(I))+real(V(lineI) .*conj(1j*lineB2.*V(lineI) ) );
+output_args=real((V(lineI)).*conj(I));
 end
 

BranchQ.m

@@ -1,4 +1,4 @@
 function [ output_args ] = BranchQ( V,I,lineI,lineB2 )
-output_args=imag((V(lineI)).*conj(I))+imag(V(lineI) .*conj(1j*lineB2.*V(lineI) ) );
+output_args=imag((V(lineI)).*conj(I));
 end
 

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('s1047.dat', '0');
+[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee14.dat', '0');
 %% 开始生成量测量
 sigma=0.03;% 标准差
 %% 电压
@@ -27,7 +27,7 @@ lineX=newwordParameter.line.lineX;
 lineB2=newwordParameter.line.lineB2;
 lineG=real(1./(lineR+1j*lineX));
 lineB=imag(1./(lineR+1j*lineX));
-cmpBranchI=BranchI( cmpV,lineI,lineJ,lineR,lineX );%复数支路电流
+cmpBranchI=BranchI( cmpV,lineI,lineJ,lineR,lineX,lineB2 );%复数支路电流
 rBranchI=abs(cmpBranchI);% 支路电流幅值
 mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量
 %% 支路功率
@@ -232,13 +232,16 @@ 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
+    % 二阶导数
     
-    SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngle),lineI,lineJ,lineR,lineX );%复数支路电流
+    
+    SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngle),lineI,lineJ,lineR,lineX,lineB2 );%复数支路电流
     SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
     SEBranchQ=BranchQ( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
     SETransP=TransPower( newwordParameter,SEVolt,SEVAngle );

公式/公式.tex

@@ -61,11 +61,25 @@ Q_{ij}&=-[V_1^2-V_1V_2cos(\theta_1 - \theta_2)]B_{ij}-V_1V_2sin (\theta_1 - \the
 &=-V_1^2B_{ij}-V_1V_2[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
 \end{aligned}
 \end{equation}
+考虑接地支路后(利用Ali Abur 书上的公式)
+\newline
+有功传输功率
+\begin{equation}
+\begin{aligned}
+P_{ij}&=V_1^2(G_{ij}+G_{is})-V_1V_2[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
+\end{aligned}
+\end{equation}
+无功传输功率
+\begin{equation}
+\begin{aligned}
+Q_{ij}&=-V_1^2(B_{ij}+B_{is})-V_1V_2[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
+\end{aligned}
+\end{equation}
 线路有功功率Jacobi
 \begin{equation}
 \begin{aligned}
 \frac{\partial P_{ij}}{\partial V_1}=
-2V_1G_{ij}-V_2[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
+2V_1(G_{ij}+G_{is})-V_2[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
 \end{aligned}
 \end{equation}
 
@@ -94,7 +108,7 @@ Q_{ij}&=-[V_1^2-V_1V_2cos(\theta_1 - \theta_2)]B_{ij}-V_1V_2sin (\theta_1 - \the
 \begin{equation}
 \begin{aligned}
 \frac{\partial Q_{12}}{\partial V_1}&=
--2V_1B_{12}-V_2[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
+-2V_1(B_{12}+B_{is})-V_2[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
 \end{aligned}
 \end{equation}
 
@@ -206,20 +220,21 @@ ps.已检验过线路的公式。
 以上公式已经可以完成状态估计，若要实现更好的收敛性，需要利用二阶导数。
 \par
 线路支路功率二阶导数
+有功部分
 \begin{equation}
 \begin{aligned}
 \frac{\partial^2 P_{12}}{\partial V_1^2}&=
-\frac{-2}{k^2}B_{12}\\
+\frac{2}{k^2}(G_{12}+G_{1s})\\
-&=\frac{-2B_{12}}{k^2}
+&=\frac{2 (G_{12}+G_{1s} ) } {k^2}
 \end{aligned}
 \end{equation}
 
 \begin{equation}
 \begin{aligned}
 \frac{\partial^2 P_{12}}{\partial V_1 \partial V_2 }&=
-\frac{-1}{k}[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}] \\
+\frac{-1}{k}[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}] \\
 &=
-\frac{-[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]}{k}
+\frac{-[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]}{k}
 \end{aligned}
 \end{equation}
 
@@ -232,26 +247,69 @@ ps.已检验过线路的公式。
 \begin{equation}
 \begin{aligned}
 \frac{\partial P_{12}}{\partial \theta_1^2}&=
-V_1V_2[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
+\frac{V_1}{k} V_2[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
 \end{aligned}
 \end{equation}
 
 \begin{equation}
 \begin{aligned}
 \frac{\partial P_{12}}{\partial \theta_1 \partial \theta_2}&=
-V_1V_2[-cos(\theta_1 - \theta_2)G_{ij}-sin (\theta_1 - \theta_2)B_{ij}] \\
+\frac{V_1}{k}  V_2[-cos(\theta_1 - \theta_2)G_{ij}-sin (\theta_1 - \theta_2)B_{ij}] \\
 &=
--V_1V_2[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
+- \frac{V_1}{k} V_2[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
 \end{aligned}
 \end{equation}
 
 \begin{equation}
 \begin{aligned}
 \frac{\partial P_{12}}{\partial \theta_2^2}&=
--V_1V_2[-cos(\theta_1 - \theta_2)G_{ij}-sin (\theta_1 - \theta_2)B_{ij}] \\
+- \frac{V_1}{k} V_2[-cos(\theta_1 - \theta_2)G_{ij}-sin (\theta_1 - \theta_2)B_{ij}] \\
 &=
-V_1V_2[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
+\frac{V_1}{k} V_2[cos(\theta_1 - \theta_2)G_{ij}+sin (\theta_1 - \theta_2)B_{ij}]
+\end{aligned}
+\end{equation}
+无功部分
+\begin{equation}
+\begin{aligned}
+\frac{\partial^2 Q_{12}}{\partial V_1 \partial V_1}&=
+-2\frac{( B_{12}+B_{1s})}{k^2}
 \end{aligned}
 \end{equation}
 
+\begin{equation}
+\begin{aligned}
+\frac{\partial^2 Q_{12}}{\partial V_1 \partial V_2}&=
+-\frac{1}{k}[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
+\end{aligned}
+\end{equation}
+
+\begin{equation}
+\begin{aligned}
+\frac{\partial^2 Q_{12}}{\partial V_2^2}&=0
+\end{aligned}
+\end{equation}
+
+\begin{equation}
+\begin{aligned}
+\frac{\partial^2 Q_{12}}{\partial \theta_1^2}&=
+-\frac{V_1}{k}V_2[-sin(\theta_1 - \theta_2)G_{ij}+cos (\theta_1 - \theta_2)B_{ij}] \\
+&=\frac{V_1}{k}V_2[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
+\end{aligned}
+\end{equation}
+
+\begin{equation}
+\begin{aligned}
+\frac{\partial^2 Q_{12}}{\partial \theta_1 \partial \theta_2}&=
+-\frac{V_1}{k}V_2[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
+\end{aligned}
+\end{equation}
+
+\begin{equation}
+\begin{aligned}
+\frac{\partial^2 Q_{12}}{\partial \theta_2^2}
+&=\frac{V_1}{k}V_2[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
+\end{aligned}
+\end{equation}
+以上公式对线路和没有计及接地支路的变压器适用，只是线路中变比$k$为1.
+
 \end{document}