From 23d116883beb84ebb910280aa6808ee281b68006 Mon Sep 17 00:00:00 2001
From: "dmy@lab" <dmy@lab.lab>
Date: Thu, 4 Jun 2015 22:38:28 +0800
Subject: [PATCH] =?UTF-8?q?=E8=AE=A1=E7=AE=97=E6=97=B6=E9=97=B4=E6=AF=94?=
 =?UTF-8?q?=E8=BE=83=E5=A4=9A=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Signed-off-by: dmy@lab <dmy@lab.lab>
---
 OPF.m        |  7 ++++---
 SolveIt.m    | 47 ++++++++++++++++++++++++++++-------------------
 jacobian_M.m |  8 ++++----
 3 files changed, 36 insertions(+), 26 deletions(-)

diff --git a/OPF.m b/OPF.m
index 7b830b5..6338a3c 100644
--- a/OPF.m
+++ b/OPF.m
@@ -309,16 +309,16 @@ while(maxD>Precision)
     end
     %% 开始计算OPF
     %% 形成等式约束的雅克比
-    dH_dx = jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);  %形成雅克比矩阵
+%     dH_dx = jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);  %形成雅克比矩阵
     %% 形成不等式约束的雅克比
     %%  开始构建ddg
     %%  开始构建deltF
 %     deltF=func_deltF(wVolt,wPD,wQD,mPD3P,PD3P,QD3P,mQD3P,Volt,mVoltABCV,Busnum,Loadi);
     %% 形成方程矩阵
-    Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi);
+%     Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi);
     %% 开始解方程
 %         fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
-    XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,mPD3P,PD3P,mQD3P,QD3P,mVoltABCV,Volt,Mat_H,wVolt,wPD,wQD,dH_dx);
+    XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,mPD3P,PD3P,mQD3P,QD3P,mVoltABCV,Volt,wVolt,wPD,wQD,Y,Yangle,UAngel,r,c);
     %%取各分量
     [deltX,deltY]=AssignXX(XX,ContrlCount,Loadi,Balance,Busnum);
     [Init_Y,PG,QG,Volt,UAngel,PD3P,QD3P]=Modification(Init_Y,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD3P,QD3P,Loadi);
@@ -328,6 +328,7 @@ while(maxD>Precision)
 end
 totalTime=toc
 (rVoltABCV-Volt)./rVoltABCV*100;
+Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi);
 PD3P=-Mat_H(Loadi);
 QD3P=-Mat_H(Loadi+Busnum*3);
 %% 计算统计量
diff --git a/SolveIt.m b/SolveIt.m
index 3fa66a8..73d6a46 100644
--- a/SolveIt.m
+++ b/SolveIt.m
@@ -1,24 +1,33 @@
-function XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,PD0,PD,QD0,QD,mVolt,Volt,Mat_H,wVolt,wPD,wQD,dH_dx)
+function XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,PD0,PD,QD0,QD,mVolt,Volt,wVolt,wPD,wQD,Y,Yangle,UAngel,r,c)
 % dH_dx
-dP_x=-dH_dx(Loadi,:);
-dQ_x=-dH_dx(Loadi+Busnum*3,:);
-dV_x=[sparse(1:Busnum*3,1:Busnum*3,1,Busnum*3,Busnum*3),sparse(Busnum*3,Busnum*3)];
+dP_x = jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);  %形成雅克比矩阵
+dP_x=-dP_x(Loadi,:);
+dQ_x = jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);
+dQ_x=-dQ_x(Loadi+Busnum*3,:);
+dV_x=[eye(Busnum*3,Busnum*3),zeros(Busnum*3,Busnum*3)];
 H=[dP_x;dQ_x;dV_x];
-c=dH_dx(setdiff(1:Busnum*3,Loadi),:);
-c=[c;dH_dx(setdiff(1:Busnum*3,Loadi)+Busnum*3,:)];
-aa=[H'*diag([wPD;wQD;wVolt])*H c'
-    c zeros(size(c,1),size(c',2))];
-% aa=H'*eye(length(H))*H;
-PD=-Mat_H(Loadi,:);
-QD=-Mat_H(Loadi+Busnum*3,:);
+cP=jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);
+cP=cP(setdiff(1:Busnum*3,Loadi),:);
+cQ=jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);
+cQ=cQ(setdiff(1:Busnum*3,Loadi)+Busnum*3,:);
+C=[cP;cQ];
+aa=[H'*diag([wPD;wQD;wVolt])*H C'
+    C zeros(size(C,1),size(C',2))];
+AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Yangle,Busnum*3,Busnum*3);
+PD=-diag(Volt)*Y.*cos(AngleIJ)*Volt;
+PD=PD(Loadi);
+QD=-diag(Volt)*Y.*sin(AngleIJ)*Volt;
+QD=QD(Loadi);
+dP=diag(Volt)*Y.*cos(AngleIJ)*Volt;
+dQ=diag(Volt)*Y.*sin(AngleIJ)*Volt;
 deltZ=[(PD0-PD);
     (QD0-QD);
     (mVolt-Volt);
     ];
 yy=[
     H'*diag([wPD;wQD;wVolt])*deltZ;
-    -Mat_H(setdiff(1:Busnum*3,Loadi));
-    -Mat_H(setdiff(1:Busnum*3,Loadi)+Busnum*3);
+    -dP(setdiff(1:Busnum*3,Loadi));
+    -dQ(setdiff(1:Busnum*3,Loadi));
 ];
 % yy=[
 %     H'*eye(length(H))*deltZ;
@@ -31,9 +40,9 @@ aa(t+(Balance-1)*3+3,:)=0;
 aa(:,t+(Balance-1)*3+1)=0;
 aa(:,t+(Balance-1)*3+2)=0;
 aa(:,t+(Balance-1)*3+3)=0;
-aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
-aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
-aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
 
 % aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
 % aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
@@ -47,9 +56,9 @@ aa(t+(Balance-1)*3+3,:)=0;
 aa(:,t+(Balance-1)*3+1)=0;
 aa(:,t+(Balance-1)*3+2)=0;
 aa(:,t+(Balance-1)*3+3)=0;
-aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
-aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
-aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
 % aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
 % aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
 % aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
diff --git a/jacobian_M.m b/jacobian_M.m
index 90cfa57..1a3c74e 100644
--- a/jacobian_M.m
+++ b/jacobian_M.m
@@ -8,10 +8,10 @@ function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c)
 %% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
 AngleIJ=UAngel(r)-UAngel(c)-Angle;
 mat_AngleIJ=sparse(r,c,AngleIJ,Busnum*3,Busnum*3);
-H=diag(Volt)*Y.*sin(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)-sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
-N=-diag(Volt)*Y.*cos(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)+sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
-J=sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*cos(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
-L=sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*sin(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
+H=diag(Volt)*Y.*sin(mat_AngleIJ)*diag(Volt)-diag(Y.*sin(mat_AngleIJ)*Volt)*diag(Volt);
+N=-diag(Volt)*Y.*cos(mat_AngleIJ)*diag(Volt)+diag(Y.*cos(mat_AngleIJ)*Volt)*diag(Volt);
+J=diag(Y.*cos(mat_AngleIJ)*Volt)+Y.*cos(mat_AngleIJ)*diag(Volt);
+L=diag(Y.*sin(mat_AngleIJ)*Volt)+Y.*sin(mat_AngleIJ)*diag(Volt);
 t1=[J,H;
     L,N;
 ];