From 3dd10280c7606b2a34cf6d1b8119272026163fe6 Mon Sep 17 00:00:00 2001
From: "dmy@lab" <dmy@lab.lab>
Date: Tue, 5 May 2015 22:35:55 +0800
Subject: [PATCH] =?UTF-8?q?=E4=BF=AE=E6=94=B9=E4=B8=BA=E5=8F=AA=E5=8F=AA?=
 =?UTF-8?q?=E6=9C=89=E4=B8=8D=E5=AD=98=E5=9C=A8=E5=AE=9E=E6=97=B6=E6=95=B0?=
 =?UTF-8?q?=E6=8D=AE=E7=9A=84=E6=97=B6=E5=80=99=E6=89=8D=E5=8A=A0=E5=88=B0?=
 =?UTF-8?q?=E4=B8=8D=E7=AD=89=E5=BC=8F=E4=B8=AD=EF=BC=8C=E6=94=B6=E6=95=9B?=
 =?UTF-8?q?=E9=80=9F=E5=BA=A6=E5=BF=AB=E5=BE=88=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>
---
 FormGstate1.m      |  6 +++---
 IPMLoop.m          | 43 +++++++++++++++++++++++++------------------
 func_deltGstate1.m | 11 +++++++----
 run.m              | 10 +++++-----
 test.m             |  2 +-
 5 files changed, 41 insertions(+), 31 deletions(-)

diff --git a/FormGstate1.m b/FormGstate1.m
index 7972cd7..3c9a919 100644
--- a/FormGstate1.m
+++ b/FormGstate1.m
@@ -1,4 +1,4 @@
-function Mat_G=FormGstate1(I1r,I1i)
+function Mat_G=FormGstate1(I1r,I1i,Loadi)
 %t1=PG(PVi);
 %GP=t1;%发电机P
 %GP=[4.5 4.5]';
@@ -10,7 +10,7 @@ function Mat_G=FormGstate1(I1r,I1i)
 % t4=sum(t3,2);%发电机Q
 %GQ=t4;
 Mat_G=[
-    I1r;
-    I1i;
+    I1r(Loadi);
+    I1i(Loadi);
     ];
 end
\ No newline at end of file
diff --git a/IPMLoop.m b/IPMLoop.m
index 16cc370..d62d9b5 100644
--- a/IPMLoop.m
+++ b/IPMLoop.m
@@ -8,12 +8,15 @@ I1r=real(guessIf1);
 I1i=imag(guessIf1);
 KK=0;
 plotGap=zeros(1,60);
-%初始化
-%状态量为 SEPD SEQD SEVmf1 SEVaf1
-RestraintCount=length(Loadi)*2;
+
 % ContrlCount=busNum*2+length(Loadi)*2;
 noBoundedLoadi=setdiff(Loadi,noLoadi);
 % noBoundedLoadi=[];
+boundedLoadi=setdiff(Loadi,noBoundedLoadi);
+%初始化
+%状态量为 SEPD SEQD SEVmf1 SEVaf1
+% RestraintCount=length(Loadi)*2;
+RestraintCount=length(boundedLoadi)*2;
 ContrlCount=busNum*2+length(Loadi)*2;
 CenterA=0.1;
 Init_Z=1*sparse(ones(RestraintCount,1));
@@ -22,7 +25,7 @@ Init_L=1*sparse(ones(RestraintCount,1));
 Init_U=1*sparse(ones(RestraintCount,1));
 Init_Y=sparse(2*busNum,1);%等式约束乘子
 Gap=(Init_L'*Init_Z-Init_U'*Init_W);
-kmax=100;
+kmax=20;
 isSetBound=0;
 largerBound=0;
 realBound=1;
@@ -37,7 +40,8 @@ while(abs(Gap)>1e-4)
     %% 形成等式约束的雅克比
     deltH=func_deltH(busNum,fsY1,Loadi,Balance);
     %% 形成不等式约束的雅克比
-     deltG=func_deltGstate1(busNum,Loadi,I1r,I1i);
+%      deltG=func_deltGstate1(busNum,Loadi,I1r,I1i);
+deltG=func_deltGstate1(busNum,Loadi,boundedLoadi,I1r,I1i);
     %%
 %     L_1Z=diag(Init_Z./Init_L);
 %     U_1W=diag(Init_W./Init_U);
@@ -55,22 +59,24 @@ while(abs(Gap)>1e-4)
     %%
     Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1);
     Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1);
-    Mat_G=FormGstate1(I1r,I1i);
+%     Mat_G=FormGstate1(I1r,I1i);
+    Mat_G=FormGstate1(I1r,I1i,ismember(Loadi,boundedLoadi) );
     Mat_H=FormH(fsY1,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance);
     Ly=Mat_H;
     
 %     if isSetBound==0
     
-        Ir=real(guessIf1);
+        Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) );
         pIr=find(Ir>0);
         nIr=find(Ir<0);
-        Ii=imag(guessIf1);
+        Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
         pIi=find(Ii>0);
         nIi=find(Ii<0);
-        lower=ones(length(Loadi)*2,1);
+%         lower=ones(length(Loadi)*2,1);
+lower=ones(length(boundedLoadi)*2,1);
 %     if abs(Gap)<1 || isSetBound==1
 
-%     if abs(Gap)<0.01
+%     if abs(Gap)<0.001
 %         largerBound=0;
 %         realBound=1;
 %     end
@@ -85,8 +91,8 @@ while(abs(Gap)>1e-4)
         lower(nIr)=(1.300)*Ir(nIr);
         lower(pIi+length(Ir))=(0.7)*Ii(pIi);
         lower(nIi+length(Ir))=(1.30)*Ii(nIi);
-        lower( ismember(Loadi,noBoundedLoadi))=-0.9;
-        lower(ismember(Loadi,noBoundedLoadi)+length(Ir))=-0.9;
+%         lower( ismember(Loadi,noBoundedLoadi))=-1.9;
+%         lower(ismember(Loadi,noBoundedLoadi)+length(Ir))=-1.9;
      end
         
     if largerBound==1
@@ -112,13 +118,13 @@ while(abs(Gap)>1e-4)
 %         I1r(I1r./lowerR>0.9998)=lowerR(I1r./lowerR>0.9998);
 % I1i(I1i./lowerI>0.9998)=lowerI(I1i./lowerI>0.9998);
     
-        Ir=real(guessIf1);
+        Ir=real(guessIf1(ismember(Loadi,boundedLoadi)));
         pIr=find(Ir>0);
         nIr=find(Ir<0);
-        Ii=imag(guessIf1);
+        Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
         pIi=find(Ii>0);
         nIi=find(Ii<0);
-        upper=ones(length(Loadi)*2,1);
+        upper=ones(length(boundedLoadi)*2,1);
 %     if abs(Gap)<1 || isSetBound==1
    if realBound==1
 %         upper(pIr)=(1.10+3/(KK+1))*Ir(pIr);
@@ -130,8 +136,8 @@ while(abs(Gap)>1e-4)
         upper(nIr)=(0.7)*Ir(nIr);
         upper(pIi+length(Ir))=(1.30)*Ii(pIi);
         upper(nIi+length(Ir))=(0.7)*Ii(nIi);
-        upper(ismember(Loadi,noBoundedLoadi))=0.9;
-        upper(ismember(Loadi,noBoundedLoadi)+length(Ir))=0.9;
+%         upper(ismember(Loadi,noBoundedLoadi))=1.9;
+%         upper(ismember(Loadi,noBoundedLoadi)+length(Ir))=1.9;
    end
 
     if largerBound==1
@@ -160,13 +166,14 @@ while(abs(Gap)>1e-4)
     [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum,Balance);
     [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi,Vref);
     Gap=(Init_L'*Init_Z-Init_U'*Init_W);+max(abs(deltX));
+%     Gap=max([max(abs(Mat_H)),max(Lz),max(Lw) ]);
     fprintf('Gap %f :%d\n',full(Gap),KK+1);
     KK=KK+1;
 end
 toc
 calcuTime=toc/KK;
 % f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(V1measurement)-V1r(Loadi)).^2)+sum((imag(V1measurement)-V1i(Loadi)).^2);
-if abs(Gap)<0.00001 && KK<kmax
+if abs(Gap)<1e-4 && KK<kmax
     isConverged=1;
 else
     isConverged=0;
diff --git a/func_deltGstate1.m b/func_deltGstate1.m
index a4aab6e..0111843 100644
--- a/func_deltGstate1.m
+++ b/func_deltGstate1.m
@@ -1,7 +1,10 @@
-function deltG=func_deltGstate1(busNum,Loadi,I1r,I1i)
+function deltG=func_deltGstate1(busNum,Loadi,boundedLoadi,I1r,I1i)
+ind_=find(ismember(Loadi,boundedLoadi));
 deltG=[
-    zeros(busNum*2,length(Loadi)*2);
-    1*eye(length(Loadi)*1),zeros(length(Loadi)*1);%I1r
-    zeros(length(Loadi)*1),1*eye(length(Loadi)*1);%I1i
+    zeros(busNum*2,length(boundedLoadi)*2);
+    sparse(ind_,1:length(boundedLoadi),1,length(I1r),length(boundedLoadi)),zeros(length(I1i)*1,length(boundedLoadi));
+    zeros(length(I1r)*1,length(boundedLoadi)),sparse(ind_,1:length(boundedLoadi),1,length(I1i),length(boundedLoadi));
+%     1*eye(length(Loadi)*1),zeros(length(I1i)*1);%I1r
+%     zeros(length(Loadi)*1),1*eye(length(Loadi)*1);%I1i
 ];
 end
\ No newline at end of file
diff --git a/run.m b/run.m
index d06df67..72bbf9e 100644
--- a/run.m
+++ b/run.m
@@ -4,10 +4,10 @@ function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_m
 close all
 clc
 clear
-lineZ=readLineZ('feeder26\lineParameter.txt');
+lineZ=readLineZ('feeder13\lineParameter.txt');
 [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
-    cap]=dataRead(lineZ,'feeder26\data1.txt');
+    cap]=dataRead(lineZ,'feeder13\data1.txt');
 a=exp(1j*2*pi/3);
 Tp2f=1/3*[1 1 1;
     1 a a^2;
@@ -291,11 +291,11 @@ wV0i=wV0i(Loadi);
 
 %% 不同量测配置
 noLoadi=[1,8,10,12];
-noLoadi=[1:7:max(Loadi)];
-noLoadi=[1,19];
+% noLoadi=[1:3:max(Loadi)];
+% noLoadi=[1,19];
 % noLoadi=[4,3,2,5,6,7,11,14,22];
 % noLoadi=Loadi;
-% noLoadi=[];
+noLoadi=[1];
 wI0r(ismember(Loadi,noLoadi))=1e100000;
 wI0i(ismember(Loadi,noLoadi))=1e100000;
 wI1r(ismember(Loadi,noLoadi))=1e100000;
diff --git a/test.m b/test.m
index c68f3a5..3f3eb2b 100644
--- a/test.m
+++ b/test.m
@@ -14,7 +14,7 @@ totalTime2Sum=10000;
 totalTime0Sum=10000;
 timeOfPrepare=0.00272333333;%准备数据的时间
 timeOfPrepare=0;
-N=500;
+N=200;
 loopN=1;
 while 1
     [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0]=run();