1.更新了公式

2.修复了dQij_dVi的公式

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

run.m

@@ -33,6 +33,8 @@ lineX=newwordParameter.line.lineX;
 lineB2=newwordParameter.line.lineB2;
 lineG=real(1./(lineR+1j*lineX));
 lineB=imag(1./(lineR+1j*lineX));
+% 处理线路电阻或电抗为0的情况，即消除NaN
+zerosRXInd=union(find(abs(lineR)<1e-5),find(abs(lineX)<1e-5));
 cmpBranchI=BranchI( cmpV,lineI,lineJ,lineR,lineX );%复数支路电流
 rBranchI=abs(cmpBranchI);% 支路电流幅值
 mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量
@@ -47,6 +49,9 @@ mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%֧·
 rBranchQ=BranchQ( cmpV,cmpBranchI,lineI,lineB2 );
 % rBranchQ=rBranchQ(abs(rBranchQ)>1e-5);
 mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量
+% 处理线路电阻或电抗为0的情况，即消除NaN
+mBranchP=mBranchP(setdiff( 1:length(mBranchP),zerosRXInd));
+mBranchQ=mBranchQ(setdiff( 1:length(mBranchQ),zerosRXInd));
 %% 注入功率
 rPD=PD;
 PDi=find(PD~=0);
@@ -106,7 +111,7 @@ Iteration=0;
 while maxD>1e-3;
 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),lineJ, ...
+dLPij_dVi=sparse(1:length(lineI),lineI, ...
     -SEVolt(lineJ).*( ...
     lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
     )... 
@@ -117,7 +122,7 @@ dLPij_dVj=sparse(1:length(lineI),lineJ, ...
     lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
     ) ...
     ,length(lineI),length(mVolt));%线路的
-dLPij_dThetai=sparse(1:length(lineI),lineJ, ...
+dLPij_dThetai=sparse(1:length(lineI),lineI, ...
     SEVolt(lineI).*SEVolt(lineJ).*( ...
     lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
     )... 
@@ -128,9 +133,9 @@ dLPij_dThetaj=sparse(1:length(lineI),lineJ, ...
     )... 
     ,length(lineI),length(mVolt));%线路的
 
-dLQij_dVi=sparse(1:length(lineI),lineJ, ...
+dLQij_dVi=sparse(1:length(lineI),lineI, ...
     -SEVolt(lineJ).*( ...
-    lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
+    lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
     )... 
     -2*(lineB+lineB2).*SEVolt(lineI) ...
     ,length(lineI),length(mVolt));%线路的
@@ -139,7 +144,7 @@ dLQij_dVj=sparse(1:length(lineI),lineJ, ...
     lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
     ) ...
     ,length(lineI),length(mVolt));%线路的
-dLQij_dThetai=sparse(1:length(lineI),lineJ, ...
+dLQij_dThetai=sparse(1:length(lineI),lineI, ...
     -SEVolt(lineI).*SEVolt(lineJ).*( ...
     lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
     ) ...
@@ -150,7 +155,6 @@ dLQij_dThetaj=sparse(1:length(lineI),lineJ, ...
     ) ...
     ,length(lineI),length(mVolt));%线路的
 % 处理线路电阻或电抗为0的情况，即消除NaN
-zerosRXInd=union(find(abs(lineR)<1e-5),find(abs(lineX)<1e-5));
 dLPij_dVi=dLPij_dVi(setdiff( 1:size(dLPij_dVi,1),zerosRXInd),:);
 dLPij_dVj=dLPij_dVj(setdiff( 1:size(dLPij_dVj,1),zerosRXInd),:);
 dLPij_dThetai=dLPij_dThetai(setdiff( 1:size(dLPij_dThetai,1),zerosRXInd),:);
@@ -160,16 +164,26 @@ dLQij_dVj=dLQij_dVj(setdiff( 1:size(dLQij_dVj,1),zerosRXInd),:);
 dLQij_dThetai=dLQij_dThetai(setdiff( 1:size(dLQij_dThetai,1),zerosRXInd),:);
 dLQij_dThetaj=dLQij_dThetaj(setdiff( 1:size(dLQij_dThetaj,1),zerosRXInd),:);
 % 对量测值做同样处理
-mBranchP=mBranchP(setdiff( 1:length(mBranchP),zerosRXInd));
+
-mBranchQ=mBranchQ(setdiff( 1:length(mBranchQ),zerosRXInd));
 %% 进入迭代
-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;
-  dLQij_dVi+dLQij_dVj, dLQij_dThetai+dLQij_dThetaj];%jacobi
+%   dLQij_dVi+dLQij_dVj, dLQij_dThetai+dLQij_dThetaj];%jacobi
-% H=[dV_dV];
+
-h=[SEVolt;mBranchP;mBranchQ];
+H=[dV_dV;
-z=[mVolt;mBranchP;mBranchQ];
+    dLPij_dVi+dLPij_dVj ;
-W=sparse(1:length(mVolt),1:length(mVolt),1/sigma.^2,length(mVolt),length(mVolt));
+    dLQij_dVi+dLQij_dVj ];%jacobi
+
+SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngel),lineI,lineJ,lineR,lineX );%复数支路电流
+SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngel),SEBranchI,lineI,lineB2 );
+SEBranchQ=BranchQ( SEVolt.*exp(1j*SEVAngel),SEBranchI,lineI,lineB2 );
+SEBranchP=SEBranchP(setdiff( 1:length(SEBranchP),zerosRXInd));
+SEBranchQ=SEBranchQ(setdiff( 1:length(SEBranchQ),zerosRXInd));
+% h=[SEVolt;SEBranchP;SEBranchQ];
+% z=[mVolt;mBranchP;mBranchQ];
+h=[SEVolt;SEBranchP;SEBranchQ];
+z=[mVolt;mBranchP;mBranchP];
+W=sparse(1:length(h),1:length(h),1/sigma.^2,length(h),length(h));
 G=H'*W*H;
 g=-H'*W*(z-h);
 % 平衡节点相角恒定;
@@ -182,11 +196,10 @@ maxD=max(abs(dX))
 % 更新变量
 SEVolt=SEVolt+dX(1:length(mVolt));
 Iteration=Iteration+1;
-SEVAngel=SEVAngel+dX(length(mVolt)+1:end);
+% SEVAngel=SEVAngel+dX(length(mVolt)+1:end);
 end
 %% 输出结果
 fprintf('迭代%d次\n',Iteration);
-h=SEVolt;
 fval=full((z-h)'*W*(z-h));
 fprintf('目标函数为:%f\n',fval);
 fprintf('相对误差\n');

公式/公式.tex

@@ -25,7 +25,7 @@
 \end{equation}
 \begin{equation}
 \begin{aligned}
-\dot{S}^*_{12}&=V_1e^{j \theta_1}
+\dot{S}_{12}&=V_1e^{j \theta_1}
 \dot{I}^*_{12} \\
 &=V_1e^{j \theta_1}
 (V_1e^{-j \theta_1} - V_2e^{-j \theta_2})
@@ -50,15 +50,74 @@
 有功传输功率
 \begin{equation}
 \begin{aligned}
-\dot{P}_{ij}&=[V_1^2-V_1V_2cos(\theta_1 - \theta_2)]G_{ij}-V_1V_2sin (\theta_1 - \theta_2)B_{ij} \\
+P_{ij}&=[V_1^2-V_1V_2cos(\theta_1 - \theta_2)]G_{ij}-V_1V_2sin (\theta_1 - \theta_2)B_{ij} \\
 &=V_1^2G_{ij}-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}
-\dot{Q}_{ij}&=-[V_1^2-V_1V_2cos(\theta_1 - \theta_2)]B_{ij}-V_1V_2sin (\theta_1 - \theta_2)G_{ij} \\
+Q_{ij}&=-[V_1^2-V_1V_2cos(\theta_1 - \theta_2)]B_{ij}-V_1V_2sin (\theta_1 - \theta_2)G_{ij} \\
 &=-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}
+线路有功功率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}]
+\end{aligned}
+\end{equation}
+
+\begin{equation}
+\begin{aligned}
+\frac{\partial P_{12}}{\partial V_2}=
+-V_1[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}&=
+-V_1V_2[-sin(\theta_1 - \theta_2)G_{ij}+cos (\theta_1 - \theta_2)B_{ij}] \\
+&=V_1V_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 P_{12}}{\partial \theta_2}&=
+-V_1V_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 Q_{12}}{\partial V_1}&=
+-2V_1B_{12}-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 Q_{12}}{\partial V_2}&=
+-V_1[sin(\theta_1 - \theta_2)G_{ij}-cos (\theta_1 - \theta_2)B_{ij}]
+\end{aligned}
+\end{equation}
+
+\begin{equation}
+\begin{aligned}
+\frac{\partial Q_{12}}{\partial \theta_1}&=
+-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}
+\frac{\partial Q_{12}}{\partial \theta_2}&=
+-V_1V_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}]
+\end{aligned}
+\end{equation}
+
 \end{document}