1.更新变压器功率的函数

2.删掉无用代码，更易阅读

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

TransPower.m

@@ -1,11 +1,13 @@
-function [ output_args ] = TransPower( newwordParameter,cmpY,Volt,VAngel )
+function [ output_args ] = TransPower( newwordParameter,Volt,VAngel )
 %TRANSPOWER Summary of this function goes here
 %   Detailed explanation goes here
 transI=newwordParameter.trans.transI;
 transJ=newwordParameter.trans.transJ;
+transR=newwordParameter.trans.transR;
+transX=newwordParameter.trans.transX;
 cmpSEV=Volt.*exp(1j*VAngel); %复数电压
 % cmpY=sparse(r,c,Y.*(1j*YAngel),length(Volt),length(Volt));
-cmpSETransI=(cmpSEV(transI)-cmpSEV(transJ)).*cmpY(sub2ind(size(cmpY),transI,transJ));
+cmpSETransI=(cmpSEV(transI)-cmpSEV(transJ)).*1./(transR+1j*transX);
 output_args=real(cmpSEV(transI).*conj(cmpSETransI));
 end
 

TransReactivePower.m

@@ -1,11 +1,13 @@
-function [ output_args ] = TransReactivePower( newwordParameter,cmpY,Volt,VAngel )
+function [ output_args ] = TransReactivePower( newwordParameter,Volt,VAngel )
 %TRANSREAVTIVEPOWER Summary of this function goes here
 %   Detailed explanation goes here
 transI=newwordParameter.trans.transI;
 transJ=newwordParameter.trans.transJ;
+transR=newwordParameter.trans.transR;
+transX=newwordParameter.trans.transX;
 cmpSEV=Volt.*exp(1j*VAngel); %复数电压
 % cmpY=sparse(r,c,Y.*(1j*YAngel),length(Volt),length(Volt));
-cmpSETransI=(cmpSEV(transI)-cmpSEV(transJ)).*cmpY(sub2ind(size(cmpY),transI,transJ));
+cmpSETransI=(cmpSEV(transI)-cmpSEV(transJ)).*1./(transR+1j*transX);
 output_args=imag(cmpSEV(transI).*conj(cmpSETransI));
 end
 

run.m

@@ -33,25 +33,14 @@ lineX=newwordParameter.line.lineX;
 lineB2=newwordParameter.line.lineB2;
 lineG=real(1./(lineR+1j*lineX));
 lineB=imag(1./(lineR+1j*lineX));
-% 处理线路电阻或电抗为0的情况，即消除NaN
-% zerosLRXInd=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);%支路电流量测量
 %% 支路功率
 rBranchP= BranchP( cmpV,cmpBranchI,lineI,lineB2 );
-% rV=abs(cmpV);
-% rA=angle(cmpV);
-% a=rV(lineI).^2.*(lineG)-rV(lineI).*rV(lineJ).*(lineG.*cos(rA(lineI)-rA(lineJ) )+lineB.*sin(rA(lineI)-rA(lineJ) ) );
-% b=-rV(lineI).^2.*(lineB+lineB2)-rV(lineI).*rV(lineJ).*(lineG.*sin(rA(lineI)-rA(lineJ) )-lineB.*cos(rA(lineI)-rA(lineJ) ) );
-% rBranchP=rBranchP(abs(rBranchP)>1e-5);
 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),zerosLRXInd));
-% mBranchQ=mBranchQ(setdiff( 1:length(mBranchQ),zerosLRXInd));
 %% 注入功率
 rPD=PD;
 PDi=find(PD~=0);
@@ -76,12 +65,8 @@ transR=newwordParameter.trans.transR;
 transX=newwordParameter.trans.transX;
 transG=real(1./(transR+1j*transX));
 transB=imag(1./(transR+1j*transX));
-% 处理线路电阻或电抗为0的情况，即消除NaN
+rTransP=TransPower( newwordParameter,rVolt,rVAngel );
-% zerosTRXInd=union(find(abs(transR)<1e-5),find(abs(transX)<1e-5));
+rTransQ=TransReactivePower( newwordParameter,rVolt,rVAngel );
-rTransP=TransPower( newwordParameter,cmpY,rVolt,rVAngel );
-rTransQ=TransReactivePower( newwordParameter,cmpY,rVolt,rVAngel );
-% rTransP=rTransP( setdiff( 1:length(rTransP),zerosTRXInd) );
-% rTransQ=rTransQ( setdiff( 1:length(rTransQ),zerosTRXInd) );
 mTransP=rTransP.*(normrnd(0,sigma,length(rTransP),1)+1);
 mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1);
 %% 发电机注入功率
@@ -210,33 +195,20 @@ while maxD>1e-3;
         transG.*cos(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))...
         ) ...
         ,length(transI),length(mVolt));%变压器
-    % 处理线路电阻或电抗为0的情况，即消除NaN
-%     dLPij_dVi=dLPij_dVi(setdiff( 1:size(dLPij_dVi,1),zerosLRXInd),:);
-%     dLPij_dVj=dLPij_dVj(setdiff( 1:size(dLPij_dVj,1),zerosLRXInd),:);
-%     dLPij_dThetai=dLPij_dThetai(setdiff( 1:size(dLPij_dThetai,1),zerosLRXInd),:);
-%     dLPij_dThetaj=dLPij_dThetaj(setdiff( 1:size(dLPij_dThetaj,1),zerosLRXInd),:);
-%     dLQij_dVi=dLQij_dVi(setdiff( 1:size(dLQij_dVi,1),zerosLRXInd),:);
-%     dLQij_dVj=dLQij_dVj(setdiff( 1:size(dLQij_dVj,1),zerosLRXInd),:);
-%     dLQij_dThetai=dLQij_dThetai(setdiff( 1:size(dLQij_dThetai,1),zerosLRXInd),:);
-%     dLQij_dThetaj=dLQij_dThetaj(setdiff( 1:size(dLQij_dThetaj,1),zerosLRXInd),:);
-    % 对量测值做同样处理
-    
     %% 进入迭代
     % H=[dV_dV ,dV_dTyta;
     %   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,dV_dTyta;
         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 ;
+        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 );%复数支路电流
     SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngel),SEBranchI,lineI,lineB2 );
     SEBranchQ=BranchQ( SEVolt.*exp(1j*SEVAngel),SEBranchI,lineI,lineB2 );
-%     SEBranchP=SEBranchP(setdiff( 1:length(SEBranchP),zerosLRXInd));
+    SETransP=TransPower( newwordParameter,cmpY,rVolt,rVAngel );
-%     SEBranchQ=SEBranchQ(setdiff( 1:length(SEBranchQ),zerosLRXInd));
-    % 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));
@@ -250,9 +222,9 @@ while maxD>1e-3;
     dX=G\-g;
     maxD=max(abs(dX))
     % 更新变量
-    % SEVolt=SEVolt+dX(1:length(mVolt));
+    SEVolt=SEVolt+dX(1:length(mVolt));
     Iteration=Iteration+1;
-    SEVAngel=SEVAngel+dX(1:end);
+    SEVAngel=SEVAngel+dX(length(mVolt)+1:end);
 end
 %% 输出结果
 fprintf('迭代%d次\n',Iteration);