From eb222764ca97777167c13b2c5588383d570474ad Mon Sep 17 00:00:00 2001
From: unknown <dmy@dmy-PC.(none)>
Date: Tue, 11 Sep 2012 19:27:02 +0800
Subject: [PATCH] =?UTF-8?q?=E6=B7=BB=E5=8A=A0=E4=BA=86QD=E6=94=B6=E6=95=9B?=
 =?UTF-8?q?=EF=BC=8C=E4=BD=86=E6=B2=A1=E6=9C=89=E6=89=BE=E5=88=B0=E6=9C=80?=
 =?UTF-8?q?=E4=BC=98=E8=A7=A3=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Signed-off-by: unknown <dmy@dmy-PC.(none)>
---
 FormLw.m       |  2 +-
 FormLz.m       |  2 +-
 Modification.m |  9 +++++----
 OPF.m          | 13 +++++++------
 OPF_Init.m     |  5 +++--
 func_ddg.m     |  2 +-
 func_ddh.m     |  2 +-
 func_deltF.m   |  1 +
 func_deltG.m   |  6 ++++++
 func_deltH.m   |  3 ++-
 pf.m           |  2 +-
 11 files changed, 29 insertions(+), 18 deletions(-)

diff --git a/FormLw.m b/FormLw.m
index e7ecda5..b9f3e81 100644
--- a/FormLw.m
+++ b/FormLw.m
@@ -2,7 +2,7 @@ function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK)
 KK=999;
 PU=GenU(:,2);%发电机有功上界
 QU=PVQU(:,1);%发电机无功上界
-VoltU=(1.2+1/exp(KK))*ones(1,Busnum);
+VoltU=(1.2)*ones(1,Busnum);
 %VoltU=10*ones(1,Busnum);
 PDU=PD0(Loadi);
 PDU(PDU>0)=1.200*PDU(PDU>0);
diff --git a/FormLz.m b/FormLz.m
index 1964285..255ef59 100644
--- a/FormLz.m
+++ b/FormLz.m
@@ -2,7 +2,7 @@ function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK)
 KK=999;
 PL=GenL(:,2);%发电机有功下界
 QL=PVQL(:,1);%发电机无功下界
-VoltL=(0.8-1/exp(KK))*ones(1,Busnum);
+VoltL=(0.8)*ones(1,Busnum);
 %VoltL=-10*ones(1,Busnum);
 PDL=PD0(Loadi);
 PDL(PDL>0)=0.800*PDL(PDL>0);
diff --git a/Modification.m b/Modification.m
index 2e4b6f0..80f15e9 100644
--- a/Modification.m
+++ b/Modification.m
@@ -1,4 +1,4 @@
-function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,Loadi)
+function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi)
 AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
 %fprintf('AlphaP %f\n',full(AlphaP));
 AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
@@ -13,9 +13,10 @@ Init_Y=Init_Y+AlphaD*deltY';
 PG(PGi)=PG(PGi)+AlphaP*deltX(1:size(PGi,1));
 %QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
 QG(PVi)=QG(PVi)+AlphaP*deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
-t=deltX(size(PVi,1)+size(PGi,1)+1:size(PVi,1)+size(PGi,1)+size(Loadi,1));
-PD(Loadi)=PD(Loadi)+AlphaP*t;
-t=deltX(size(PVi,1)+size(PGi,1)+size(Loadi,1)+1:ContrlCount)';
+t=deltX(size(PVi,1)+size(PGi,1)+1:size(PVi,1)+size(PGi,1)+size(Loadi,1)*2);
+PD(Loadi)=PD(Loadi)+AlphaP*t(1:length(Loadi));
+QD(Loadi)=QD(Loadi)+AlphaP*t(length(Loadi)+1:2*length(Loadi));
+t=deltX(size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+1:ContrlCount)';
 t(Busnum+Balance)=0;
 %Volt=Volt+AlphaP*t(2:2:2*Busnum);暂时改一下20111227
 %UAngel=UAngel+AlphaP*t(1:2:2*Busnum);暂时改一下20111227
diff --git a/OPF.m b/OPF.m
index 8f405d6..b5a7bf7 100644
--- a/OPF.m
+++ b/OPF.m
@@ -2,10 +2,11 @@ tic
 clc
 clear
 [kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL,Liner,Linex,Lineb,Transforr,Transforx,Transfork0]= ...
-pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); 
+pf('ieee4PG.dat');
+%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); 
 %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); 
 %pf('c:/file31.txt'); 
-%pf('ieee118PG.dat');
+
 
 %% 计算功率因数
 %atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0));
@@ -28,7 +29,7 @@ PG0(Balance)=PGBal(Balance);
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 plotGap=zeros(1,50);
-ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
+ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+Busnum*2;
 kmax=60;
 %% 20120523 临时
 QD_NON_ZERO=QD(PD==0 & QD~=0);
@@ -72,7 +73,7 @@ while(abs(Gap)>Precision)
     XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
     %%取各分量
     [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
-    [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,Loadi);
+    [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi);
     Gap=(Init_L*Init_Z'-Init_U*Init_W');
     KK=KK+1;
 end
@@ -89,8 +90,8 @@ maxPDError=max(absPDLoad);
 disp('index');
 Loadi(absPDLoad==maxPDError);
 %% 计算总线损
-totalLoss=(PGBal(Balance)-sum(PD(Loadi)))*100;
-fprintf('总的损耗为%f(MW 有名值)\n',totalLoss);
+totalLoss=(sum(PG)-sum(PD(Loadi)))*100;
+fprintf('总的损耗为%f(MW 有名值)\n',full(totalLoss));
 %% 计算各线损
 Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel);
 toc  
diff --git a/OPF_Init.m b/OPF_Init.m
index d24a221..b8a966f 100644
--- a/OPF_Init.m
+++ b/OPF_Init.m
@@ -18,8 +18,9 @@ tQL=sparse(PVQL(:,1));% 
 PG(PGi)=(tPU+tPL)/2;
 QG(PVi)=(tQU+tQL)/2;
 wG=1*ones(size(PGi,1),1);
-randInt=randperm(size(Loadi,1));
-randPDind=randInt(1:10);
+%randInt=randperm(size(Loadi,1));
+%randPDind=randInt(1:10);
+randPDind=0;
 wD=1*ones(size(Loadi,1),1);
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
diff --git a/func_ddg.m b/func_ddg.m
index 1dbc89d..5f9d854 100644
--- a/func_ddg.m
+++ b/func_ddg.m
@@ -1,5 +1,5 @@
 function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi)
 
-ddg=sparse(size(PVi,1)+size(PGi,1)+size(Loadi,1)+1*Busnum,RestraintCount);
+ddg=sparse(size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+1*Busnum,RestraintCount);
 
 end
\ No newline at end of file
diff --git a/func_ddh.m b/func_ddh.m
index e678761..8403554 100644
--- a/func_ddh.m
+++ b/func_ddh.m
@@ -50,7 +50,7 @@ t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
 ];
 sizePGi=size(PGi,1);
 sizePVi=size(PVi,1);
-sizeLoadi=size(Loadi,1);
+sizeLoadi=size(Loadi,1)*2;
 ddh=[
     sparse(sizePGi+sizePVi+sizeLoadi,ContrlCount);
     sparse(2*Busnum,sizePVi+sizePGi+sizeLoadi),-t;
diff --git a/func_deltF.m b/func_deltF.m
index 3217b2f..8c71c7b 100644
--- a/func_deltF.m
+++ b/func_deltF.m
@@ -6,6 +6,7 @@ deltF=[
     sparse(t1);
     sparse(size(PVi,1),1);
     sparse(t2);
+    sparse(length(Loadi),1);
     sparse(2*Busnum,1);
 ];
 
diff --git a/func_deltG.m b/func_deltG.m
index 26d7ca0..637130b 100644
--- a/func_deltG.m
+++ b/func_deltG.m
@@ -19,6 +19,11 @@ dg2_dPD=sparse(size(Loadi,1),size(PVi,1));
 dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
 dg4_dPD=sparse(size(Loadi,1),Busnum);
 %% 
+dg1_dQD=sparse(size(Loadi,1),size(PGi,1));
+dg2_dQD=sparse(size(Loadi,1),size(PVi,1));
+dg3_dQD=sparse(length(Loadi),length(Loadi));
+dg4_dQD=sparse(size(Loadi,1),Busnum);
+%% 
 dg1_dx=sparse(2*Busnum,sizePGi);
 dg2_dx=sparse(2*Busnum,sizePVi);
 dg3_dx=sparse(2*Busnum,sizeLoadi);
@@ -29,5 +34,6 @@ dg4_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
 deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg4_dPg;
        dg1_dQr,dg2_dQr,dg3_dQr,dg4_dQr;
        dg1_dPD,dg2_dPD,dg3_dPD,dg4_dPD;
+       dg1_dQD,dg2_dQD,dg3_dQD,dg4_dQD;
        dg1_dx,dg2_dx,dg3_dx,dg4_dx;
 ];
\ No newline at end of file
diff --git a/func_deltH.m b/func_deltH.m
index 813c8e4..49ca864 100644
--- a/func_deltH.m
+++ b/func_deltH.m
@@ -3,6 +3,7 @@ function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi)
 dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum);
 dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,1),2*Busnum);
 dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)];
+dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)];
 dH_dx = jacobian_M(Busnum,Volt,Y,Angle,UAngel,r,c);  %形成雅克比矩阵
-deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx'];
+deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dQD;dH_dx'];
 end
\ No newline at end of file
diff --git a/pf.m b/pf.m
index 52a09fb..b585f71 100644
--- a/pf.m
+++ b/pf.m
@@ -10,7 +10,7 @@ function [kmax,Precision,Uangle,U,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Li
 tic;
 %% 读取数据文件
 [Busnum,Balance,PQstandard,Precision,Linei,Linej,Liner,Linex,Lineb,kmax,Transfori ,...
-    Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL]= openfile2(FileName);
+    Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL]= openfile(FileName);
 %% 形成节点导纳矩阵
 [GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,...
     Transforx,Transfork0,Branchi,Branchb);