From 361909bcde2a1dc9108adcfc40447a6d979248a7 Mon Sep 17 00:00:00 2001
From: "dmy@lab" <dmy@lab.lab>
Date: Wed, 1 Apr 2015 10:06:11 +0800
Subject: [PATCH] =?UTF-8?q?=E5=A6=82=E6=9E=9C=E7=94=A8=E9=87=8F=E6=B5=8B?=
 =?UTF-8?q?=E5=80=BC=E7=A1=AE=E5=AE=9A=E4=B8=8D=E7=AD=89=E5=BC=8F=E4=B8=8A?=
 =?UTF-8?q?=E4=B8=8B=E7=95=8C=E5=AE=B9=E6=98=93=E5=87=BA=E7=8E=B0=E4=B8=8D?=
 =?UTF-8?q?=E6=94=B6=E6=95=9B=E7=9A=84=E9=97=AE=E9=A2=98=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>
---
 FormLwstate1.m     | 12 ++++++------
 FormLzstate1.m     | 12 ++++++------
 IPMLoop.m          | 35 +++++++----------------------------
 feeder13/data1.txt |  6 +++---
 run.m              | 14 +++++++++++++-
 5 files changed, 35 insertions(+), 44 deletions(-)

diff --git a/FormLwstate1.m b/FormLwstate1.m
index 3274aec..2a3b788 100644
--- a/FormLwstate1.m
+++ b/FormLwstate1.m
@@ -8,18 +8,18 @@ pIi=find(Ii>0);
 nIi=find(Ii<0);
 % lower=-0.2*sparse(ones(length(Loadi)*2,1));
 upper=ones(length(Loadi)*2,1);
-upper(pIr)=1.2*Ir(pIr);
-upper(nIr)=0.8*Ir(nIr);
-upper(pIi+length(Ir))=1.2*Ii(pIi);
-upper(nIi+length(Ir))=0.8*Ii(nIi);
+upper(pIr)=1.4*Ir(pIr);
+upper(nIr)=0.6*Ir(nIr);
+upper(pIi+length(Ir))=1.4*Ii(pIi);
+upper(nIi+length(Ir))=0.6*Ii(nIi);
 
 %太小的数都要放宽一些
 tooSmall=find(abs(Ir)<0.0005);
 % upper(tooSmall)=0.99*abs(Ir(tooSmall));
-upper(tooSmall)=0.002;
+upper(tooSmall)=0.2;
 tooSmall=find(abs(Ii)<0.0005);
 % upper(tooSmall+length(Ir))=0.99*abs(Ii(tooSmall));
-upper(tooSmall+length(Ir))=0.002;
+upper(tooSmall+length(Ir))=0.2;
 
 % upper([4,5,6])=[1;1;1];
 % upper([4,5,6])=0.2*Ii(nIi);
diff --git a/FormLzstate1.m b/FormLzstate1.m
index 6a2e11e..d8907b1 100644
--- a/FormLzstate1.m
+++ b/FormLzstate1.m
@@ -7,18 +7,18 @@ pIi=find(Ii>0);
 nIi=find(Ii<0);
 % lower=-0.2*sparse(ones(length(Loadi)*2,1));
 lower=ones(length(Loadi)*2,1);
-lower(pIr)=0.8*Ir(pIr);
-lower(nIr)=1.2*Ir(nIr);
-lower(pIi+length(Ir))=0.8*Ii(pIi);
-lower(nIi+length(Ir))=1.2*Ii(nIi);
+lower(pIr)=0.6*Ir(pIr);
+lower(nIr)=1.4*Ir(nIr);
+lower(pIi+length(Ir))=0.6*Ii(pIi);
+lower(nIi+length(Ir))=1.4*Ii(nIi);
 
 %太小的数都要放宽一些
 tooSmall=find(abs(Ir)<0.0005);
 % lower(tooSmall)=-0.99*abs(Ir(tooSmall));
-lower(tooSmall)=-0.002;
+lower(tooSmall)=-0.2;
 tooSmall=find(abs(Ii)<0.0005);
 % lower(tooSmall+length(Ir))=-.99*abs(Ii(tooSmall));
-lower(tooSmall+length(Ir))=-0.002;
+lower(tooSmall+length(Ir))=-0.2;
 
 % lower=-ones(length(Ir)*2,1);
 % lower([4,5,6])=[-1;-1;-1];
diff --git a/IPMLoop.m b/IPMLoop.m
index 84fae39..3f6da55 100644
--- a/IPMLoop.m
+++ b/IPMLoop.m
@@ -1,6 +1,6 @@
 function [ V1r,V1i,I1r,I1i,isConverged ] = IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,Vref )
 %把每个序的循环写在这个函数中。其实也就是内点法循环。
-V1r=1*ones(busNum,1);
+V1r=Vref*ones(busNum,1);
 V1i=0*ones(busNum,1);
 I1r=0.1*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
 I1i=0.1*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
@@ -9,11 +9,7 @@ plotGap=zeros(1,60);
 %初始化
 %状态量为 SEPD SEQD SEVmf1 SEVaf1
 state=1;%state1表示用l<=I1r<=u 这样的约束形式
-if state==1
-    RestraintCount=length(Loadi)*2;
-else
-    RestraintCount=length(Loadi)*1;
-end
+RestraintCount=length(Loadi)*2;
 ContrlCount=busNum*2+length(Loadi)*2;
 CenterA=0.1;
 Init_Z=sparse(ones(RestraintCount,1));
@@ -33,11 +29,7 @@ while(abs(Gap)>0.00001)
     %% 形成等式约束的雅克比
     deltH=func_deltH(busNum,fsY1,Loadi,Balance);
     %% 形成不等式约束的雅克比
-    if state==1
-        deltG=func_deltGstate1(busNum,Loadi,I1r,I1i);
-    else
-        deltG=func_deltG(busNum,Loadi,I1r,I1i);
-    end
+     deltG=func_deltGstate1(busNum,Loadi,I1r,I1i);
     %%
 %     L_1Z=diag(Init_Z./Init_L);
 %     U_1W=diag(Init_W./Init_U);
@@ -48,31 +40,18 @@ while(abs(Gap)>0.00001)
 %     ddh=func_ddh(busNum,Loadi,Init_Z,Init_W);
     ddh=0;
     %%  开始构建ddg
-    if state==1
-        ddg=0;
-    else
-        ddg=func_ddg(busNum,Loadi,Init_Z,Init_W);
-    end
+    ddg=0;
     %%  开始构建deltF
     deltF=func_deltF(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,busNum,fsY1,Loadi,V1r,V1i,I1r,I1i);
 %     deltF=0;
     %%
     Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1);
     Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1);
-    if state==1
-        Mat_G=FormGstate1(I1r,I1i);
-    else
-        Mat_G=FormG(I1r,I1i);
-    end
+    Mat_G=FormGstate1(I1r,I1i);
     Mat_H=FormH(fsY1,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance);
     Ly=Mat_H;
-    if state==1
-        Lz=FormLzstate1(Loadi,Mat_G,Init_L,I1measurement);
-        Lw=FormLwstate1(Loadi,Mat_G,Init_U,I1measurement);
-    else
-        Lz=FormLz(Loadi,Mat_G,Init_L,I1measurement);
-        Lw=FormLw(Loadi,Mat_G,Init_U,I1measurement);
-    end
+    Lz=FormLzstate1(Loadi,Mat_G,Init_L,I1measurement);
+    Lw=FormLwstate1(Loadi,Mat_G,Init_U,I1measurement);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程
diff --git a/feeder13/data1.txt b/feeder13/data1.txt
index 3a1a90a..208b1fd 100644
--- a/feeder13/data1.txt
+++ b/feeder13/data1.txt
@@ -6,9 +6,9 @@
 0
 632 0 0 0
 0
-630 90 70 80 70 90 70
-631 60 30 60 30 70 30
-632 60 30 70 30 60 30
+630 80 70 60 50 40 30
+631 40 10 11 50 30 10
+632 120 90 120 80 90 70
 0
 
 
diff --git a/run.m b/run.m
index d746994..28c56c3 100644
--- a/run.m
+++ b/run.m
@@ -131,7 +131,7 @@ fprintf('
 rIf0=If0;
 rIf1=If1;
 rIf2=If2;
-sigma=0.03;
+sigma=0.01;
 iterPhaseASpotLoadP=phaseASpotLoadP;
 iterPhaseBSpotLoadP=phaseBSpotLoadP;
 iterPhaseCSpotLoadP=phaseCSpotLoadP;
@@ -206,6 +206,10 @@ wV1i=abs(imag(V1measurement)).*sigma;
 wV1i(abs(wV1i)<1e-5)=1e10;
 wI1r=abs(real(I1measurement)).*sigma;
 wI1i=abs(imag(I1measurement)).*sigma;
+
+wV1i=ones(length(wV1i),1);
+wV1r=ones(length(wV1i),1);
+
 % [ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 );
 % f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(rV1)-V1r).^2)+sum((imag(rV1)-V1i).^2);
 % fprintf('目标值 %f\n',full(f));
@@ -234,6 +238,10 @@ wV2r=abs(real(V2measurement)).*sigma;
 wV2i=abs(imag(V2measurement)).*sigma;
 wI2r=abs(real(I2measurement)).*sigma;
 wI2i=abs(imag(I2measurement)).*sigma;
+
+wV2i=ones(length(wV1i),1);
+wV2r=ones(length(wV1i),1);
+
 % [ V2r,V2i,I2r,I2i ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 );
 % f=sum(([real(I2measurement);imag(I2measurement)]-[I2r;I2i]).^2)+sum((real(rV2)-V2r).^2)+sum((imag(rV2)-V2i).^2);
 % fprintf('目标值 %f\n',full(f));
@@ -263,6 +271,10 @@ wV0i=abs(imag(V0measurement)).*sigma;
 wI0r=abs(real(I0measurement)).*sigma;
 wI0r(abs(wI0r)<1e-5)=1e10;
 wI0i=abs(imag(I0measurement)).*sigma;
+
+wV0i=ones(length(wV1i),1);
+wV0r=ones(length(wV1i),1);
+
 % matlabpool local 3
 tic
 for II=1:3