From ca24c4f8d5ce72a39899349e957faf6b360804b1 Mon Sep 17 00:00:00 2001
From: facat <dmy@dmy-PC.(none)>
Date: Thu, 24 May 2012 10:07:18 +0800
Subject: [PATCH] =?UTF-8?q?=E7=94=A8=E5=8A=9F=E7=8E=87=E5=9B=A0=E6=95=B0?=
 =?UTF-8?q?=E8=AE=A1=E7=AE=97QD=EF=BC=8C=E6=94=B6=E6=95=9B=E8=89=AF?=
 =?UTF-8?q?=E5=A5=BD=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Signed-off-by: facat <dmy@dmy-PC.(none)>
---
 1047glys.txt | 263 +++++++++++++++++++++++++++++++++++++++++++++++++++
 300glys.txt  | 175 ++++++++++++++++++++++++++++++++++
 FormH.asv    |   9 +-
 FormH.m      |   7 +-
 OPF.asv      |  46 +++++----
 OPF.m        |   6 +-
 OPF_Init.m   |   1 +
 7 files changed, 480 insertions(+), 27 deletions(-)
 create mode 100644 1047glys.txt
 create mode 100644 300glys.txt

diff --git a/1047glys.txt b/1047glys.txt
new file mode 100644
index 0000000..3bd226b
--- /dev/null
+++ b/1047glys.txt
@@ -0,0 +1,263 @@
+1.33970565959900
+1.29109217433127
+1.19028994968253
+1.19028994968253
+0.896055384571344
+1.10714871779409
+1.27395165341849
+1.27395165341849
+0.960070362405688
+0.960070362405688
+0.982793723247329
+0.982793723247329
+0.982793723247329
+0.982793723247329
+1.03037682652431
+0.950546840812075
+0.950546840812075
+0.927295218001612
+0.927295218001612
+0.982793723247329
+0.960070362405688
+0.960070362405688
+0.973619668702219
+0.973619668702219
+0.956133374872731
+0.956133374872731
+0.950546840812075
+1.40006111531961
+1.38184262319520
+1.06188763560293
+1.15310689954497
+1.24796702416514
+1.19028994968253
+1.19028994968253
+1.10714871779409
+1.15257199721567
+1.25050990019539
+1.19066238914390
+1.10029280170729
+1.22094072493780
+1.03037682652431
+1.26241071170120
+1.46591938806466
+1.17269176660352
+1.28995860217369
+1.14712740491738
+1.15710711351603
+1.29778762370819
+1.13345843504701
+-1.37340076694502
+0.912907721612687
+1.10714871779409
+0.950546840812075
+0.950546840812075
+1.10714871779409
+1.10714871779409
+1.10714871779409
+1.10714871779409
+1.10714871779409
+1.10714871779409
+1.10714871779409
+0.982793723247329
+0.982793723247329
+1.09127703480230
+1.09127703480230
+1.03037682652431
+1.03037682652431
+1.07898715080620
+1.10714871779409
+1.21202565652432
+1.32581766366803
+1.35212738092095
+1.33970565959900
+1.26987609243447
+1.30454427764397
+1.17600520709514
+1.10714871779409
+1.13838855122436
+1.19028994968253
+1.10714871779409
+1.14416883366802
+1.14660191889113
+1.24904577239825
+1.16590454050981
+1.17227388112848
+1.10714871779409
+1.14416883366802
+1.10714871779409
+1.10714871779409
+1.14416883366802
+1.10714871779409
+1.10714871779409
+1.21202565652432
+1.13838855122436
+1.24904577239825
+1.10714871779409
+1.21202565652432
+1.13838855122436
+1.16590454050981
+1.15257199721567
+1.24904577239825
+1.21202565652432
+1.22202532321099
+1.14660191889113
+1.16590454050981
+1.10714871779409
+1.19028994968253
+1.24904577239825
+1.10714871779409
+1.10714871779409
+1.24904577239825
+1.21202565652432
+1.10714871779409
+1.10714871779409
+1.24904577239825
+1.17600520709514
+1.19028994968253
+1.03037682652431
+1.10714871779409
+1.19028994968253
+1.50673316779942
+1.17737387614222
+1.47673509216691
+1.52431819544389
+1.10714871779409
+1.10714871779409
+1.45945744246823
+1.43751299806830
+1.56382779848025
+-1.48369333516315
+-1.50444867018198
+1.34731972565426
+1.40394116639511
+1.40394116639511
+0.891382521387048
+1.41042088281895
+1.19028994968253
+1.27933953231703
+1.31347261182381
+1.22318946703899
+1.18530245986557
+1.16590454050981
+0.994421106203713
+1.20959257130122
+1.40196442478378
+1.16956752779005
+1.32581766366803
+1.43824479449822
+1.14683905429100
+1.24904577239825
+0.920253080893300
+1.14751255727424
+1.27027810246307
+1.19028994968253
+1.22075143854975
+1.25270875967849
+1.19780860499484
+1.16829688044339
+1.19028994968253
+1.31019393504756
+1.15530423085401
+1.12013500071426
+0
+0
+1.10714871779409
+1.19595561143930
+1.24904577239825
+1.23629600792211
+1.55359900515355
+0.551654982528547
+-1.35264902928261
+-1.47750070505462
+-1.55986780014151
+1.56783775554554
+1.53458028472249
+1.52136665855276
+-1.53064574945385
+-1.53452622064889
+1.48253068881009
+-1.55772517584471
+-1.56150866816053
+-1.29249666778979
+1.46821647437660
+1.54845367659476
+-1.45401173132527
+-1.34074192865812
+-1.37340076694502
+1.44435582953651
+1.31712223065625
+1.28474488507758
+1.00148313569423
+1.08631839775787
+1.06369782240256
+1.13416916698136
+1.15257199721567
+1.22202532321099
+1.01504060568258
+1.29469930739946
+0
+1.12692015240925
+1.52205447548497
+1.15257199721567
+1.21202565652432
+-1.35212738092095
+-1.51585168468834
+-1.37340076694502
+0
+1.40411096215343
+-1.50695625567519
+-1.48790067446930
+-1.53747533091665
+1.56421747434188
+-1.34399747874101
+-1.47713180509594
+-1.35212738092095
+1.22524074621318
+1.19028994968253
+1.35673564323108
+1.47112767430373
+1.37340076694502
+1.11816755556904
+1.16590454050981
+1.16590454050981
+1.23755202568322
+1.48765509490646
+1.27933953231703
+1.21202565652432
+1.10714871779409
+-1.35970299357215
+-1.49556171333028
+-1.51406713708914
+-1.52340101694294
+-1.52705320464371
+-1.51277025340232
+-1.35147185823543
+-1.37340076694502
+0
+1.19028994968253
+1.52915374769631
+1.15710711351603
+1.07144960511477
+1.10714871779409
+1.20473981094748
+1.14416883366802
+1.10714871779409
+1.24904577239825
+0.689800044733608
+-1.40873215532908
+1.38932862340040
+-1.36415687581762
+-1.53198088071343
+1.26491745539004
+-1.50857146559402
+-1.42469761624386
+-0.708626272127670
+1.15257199721567
+0.376490334446523
+0.785398163397448
+1.23709415025735
+1.24904577239825
+1.07248053258459
+1.23817663522098
+-1.40564764938027
diff --git a/300glys.txt b/300glys.txt
new file mode 100644
index 0000000..ee8f72b
--- /dev/null
+++ b/300glys.txt
@@ -0,0 +1,175 @@
+1.23291813855462
+1.20658322161057
+1.05033787734375
+0.969614214282690
+1.08754338693689
+0.934288111006946
+0.844529150229879
+-1.36643570544773
+-0.188221505304771
+-0.551654982528547
+-0.512210280954117
+-1.52225150040607
+-1.33506914456117
+-0.940520438838179
+1.39408747072486
+1.44533268525224
+1.23898098497191
+1.55052883208268
+-0.977333716837938
+1.28780248391596
+-1.24153476662545
+1.23732755879586
+1.23577221571408
+1.22918083614709
+1.22918083614709
+1.23605948947808
+1.23829349860715
+1.23955238671760
+1.23555911423838
+1.23465827600933
+1.25201900446950
+1.23939976933058
+1.23861741348096
+1.23721200529510
+1.23797611386815
+1.25656442771056
+0.916713602380536
+1.24024834540447
+1.24024834540447
+1.25299832098500
+1.23899585950522
+1.24156461521103
+1.14967625241387
+1.35141207346069
+1.40006111531961
+1.21202565652432
+1.19706950682934
+1.29412798771300
+1.31019393504756
+1.38457348638376
+1.14046971367234
+1.52734543140337
+1.49948886200961
+1.11821837632420
+1.22777238637419
+1.32581766366803
+1.32971519404381
+1.40868959941175
+1.08389709498363
+1.49948886200961
+0.982793723247329
+1.16590454050981
+1.29560707766573
+0
+1.26190421578122
+1.16851261966004
+1.16590454050981
+0.953821044664279
+1.21663202730047
+1.47112767430373
+1.45004167672759
+1.18482395118491
+1.18482395118491
+1.46283330655614
+1.48103815260495
+1.24049897196564
+1.24174663195456
+1.24117191937415
+1.24904577239825
+0.976229266118710
+1.24117191937415
+1.03288728150139
+1.18224992192327
+1.07407254991928
+1.26182248128073
+1.04249077137011
+1.21202565652432
+1.12103671372506
+1.07222842115668
+1.21793722425015
+1.33394756584798
+1.25979804618936
+-1.14253217391824
+1.37008732645341
+1.13684295740100
+-0.785398163397448
+1.47112767430373
+1.11955942808299
+1.16590454050981
+1.09567694572032
+1.46835247328777
+1.55284953590561
+1.35945592748700
+1.41151231491861
+1.25231017864814
+1.43170038531283
+1.19028994968253
+1.16164017403482
+1.19302966885508
+1.32581766366803
+-1.37651699244531
+1.10714871779409
+0.896055384571344
+1.19729918068842
+1.29249666778979
+1.25603866540359
+1.15257199721567
+1.24904577239825
+0
+1.40564764938027
+1.29249666778979
+0.942000040379464
+1.06550613869550
+0.938378508000418
+1.05102793509291
+1.47354312854333
+0.984866259139967
+1.12395385830913
+1.04199070618589
+1.16590454050981
+1.12103671372506
+1.25347332094928
+1.30908282439369
+1.37178525946514
+1.55781004387472
+1.24904577239825
+1.46213511095431
+1.45518514756301
+1.34399747874101
+1.29537022309886
+1.14416883366802
+1.36947921842026
+1.25422707318249
+1.24904577239825
+1.28700221758657
+0.0216889061692343
+1.30084542593131
+1.51933590773475
+1.47393650956969
+1.31362314621978
+1.13013317523554
+1.18497641785263
+1.33641158257326
+1.27933953231703
+0.927295218001612
+1.23334025756712
+1.16759372809238
+1.11372757024711
+1.11372757024711
+1.25603866540359
+1.32298481005619
+1.42316045427414
+1.26384953610358
+1.19318290207484
+1.04272187836854
+1.39122066265548
+1.25374057358575
+0.978452216516076
+1.20906708527496
+1.48294017056526
+1.22102790509537
+1.28315368548254
+1.21073699744384
+1.49800847819097
+1.23691847280296
diff --git a/FormH.asv b/FormH.asv
index ee51ec9..5983fcf 100644
--- a/FormH.asv
+++ b/FormH.asv
@@ -12,6 +12,13 @@ t9=QG-QD;
 Mat_H(1:2:2*Busnum)=t8(1:Busnum)+t4(1:Busnum);
 Mat_H(2:2:2*Busnum)=t9(1:Busnum)+t7(1:Busnum);
 Mat_H=Mat_H';
+%% 
+QDcos=textread('300glys.txt');
+t=QD(PD==0 &&);
+aa=QD;
+QD(QD~=0)=PD(QD~=0)./tan(QDcos);
+QD(PD==0)=t;
+%% 
 %%%%一下是学姐给的公式
 AngleIJ=AngleIJMat-angle(GB);
 %dP=PG-PD-diag(Volt)*Y*cos(AngleIJ)*Volt';
@@ -19,6 +26,6 @@ dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
 dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';
 %Mat_H(1:2:2*Busnum)=dP(1:Busnum);暂时改一下 20111227
 %Mat_H(2:2:2*Busnum)=dQ(1:Busnum);暂时改一下 20111227
-Mat_H=[]'
+Mat_H=[dP;dQ;];
 
 end
\ No newline at end of file
diff --git a/FormH.m b/FormH.m
index e87423d..b43f9b5 100644
--- a/FormH.m
+++ b/FormH.m
@@ -1,4 +1,4 @@
-function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y)
+function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,QD_NON_ZERO,QD_NON_ZERO_IND)
 t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
 t2=Volt'*Volt;
 t3=t1.*t2;
@@ -12,6 +12,11 @@ t9=QG-QD;
 Mat_H(1:2:2*Busnum)=t8(1:Busnum)+t4(1:Busnum);
 Mat_H(2:2:2*Busnum)=t9(1:Busnum)+t7(1:Busnum);
 Mat_H=Mat_H';
+%% 
+QDcos=textread('1047glys.txt');
+QD(QD~=0)=PD(QD~=0)./tan(QDcos);
+QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
+%% 
 %%%%一下是学姐给的公式
 AngleIJ=AngleIJMat-angle(GB);
 %dP=PG-PD-diag(Volt)*Y*cos(AngleIJ)*Volt';
diff --git a/OPF.asv b/OPF.asv
index 36d3a2a..59e3631 100644
--- a/OPF.asv
+++ b/OPF.asv
@@ -1,21 +1,27 @@
 tic
 clear
 %[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,LineLimti,LineLimtj,LinePLimt,PG,QG,PD,QD,CenterA,LineCount,PGi,PVQU,PVQL]=pf('5sj.txt');
-[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]=pf('ieee14.dat');
+[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]=pf('ieee3001PG.dat');
 GB=full(GB);
 %PVi电压节点序号
 %PVu电压节点电压标幺值
-Volt
-UAngel*180/3.1415926
+Volt;
+UAngel*180/3.1415926;
 %sprintf('%f\n',Volt);
 %sprintf('%f\n',Angel);
 %% 初值
-
-[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL);
+PG0=PG;
+PD0=PD;
+[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0);
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 plotGap=zeros(1,50);
-ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*2;
+ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3;
+kmax=600;
+%% 20120523 临时
+QD_NON_ZERO=QD(PD==0 & QD~=0);
+QD_NON_ZERO_IND=find(PD==0 & QD~=0);
+%%
 while(abs(Gap)>Precision)
     if KK>kmax
         break;
@@ -23,8 +29,6 @@ while(abs(Gap)>Precision)
     plotGap(KK+1)=Gap;
     Init_u=Gap/2/RestraintCount*CenterA;
     AngleIJMat=repmat(UAngel',1,Busnum)-repmat(UAngel,Busnum,1);
-    %indexi=[Linei',Transfori'];
-    %indexj=[Linej',Transforj'];
     %% 开始计算OPF
     %% 形成等式约束的雅克比
     deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi);
@@ -34,42 +38,36 @@ while(abs(Gap)>Precision)
     L_1Z=diag(Init_Z./Init_L);
     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
-    deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi);
+    deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0);
     %% 形成ddHy
     ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y);
     %%  开始构建ddg
     ddg=func_ddg(PGi,PVi,Busnum,RestraintCount);
     %%  开始构建deltF
-    deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi);
+    deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum);
     
     %% 形成方程矩阵
-%    Hcoma=-deltdeltF+ddh+ddg-deltG*(L_1Z-U_1W)*deltG';
-    %AA=FormAA1(deltG,deltdeltF,ddh,ddg,deltH,Init_L,Init_U,Init_W,Init_Z,Busnum,PVi,PGi,RestraintCount,Balance);
-    %AA=FormAA(L_1Z,deltG,U_1W,Hcoma,deltH);
-    %%
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
     Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
-    Mat_G=FormG(Volt,PVi,PGi,PG,QG);
-    Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y);
+    Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD);
+    Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,QD_NON_ZERO,QD_NON_ZERO_IND);
     Ly=Mat_H;
-    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU);
+    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0);
+    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
-    %LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw);
     YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx);
-    %YY=FormYY(Init_L,Lul,Lz,Ly,Init_U,Luu,Lw,LxComa);
     %% 开始解方程
-    %XX=AA\YY;
     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);
     %%取各分量
     [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX1(XX,ContrlCount,RestraintCount,Busnum);
-    %[deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX);
-    [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel]=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);
+    [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);
     Gap=(Init_L*Init_Z'-Init_U*Init_W');
     KK=KK+1;
 end
 fprintf('迭代次数%d\n',KK);
-CalCost(GenC,PG,PGi);
+%CalCost(GenC,PG,PGi);
+ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD)
 DrawGap(plotGap);
+Volt=Volt';
 toc  
 
diff --git a/OPF.m b/OPF.m
index cfdf2a9..e9f60e6 100644
--- a/OPF.m
+++ b/OPF.m
@@ -18,6 +18,10 @@ KK=0;
 plotGap=zeros(1,50);
 ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3;
 kmax=600;
+%% 20120523 临时
+QD_NON_ZERO=QD(PD==0 & QD~=0);
+QD_NON_ZERO_IND=find(PD==0 & QD~=0);
+%%
 while(abs(Gap)>Precision)
     if KK>kmax
         break;
@@ -46,7 +50,7 @@ while(abs(Gap)>Precision)
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
     Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
     Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD);
-    Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y);
+    Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,QD_NON_ZERO,QD_NON_ZERO_IND);
     Ly=Mat_H;
     Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0);
     Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0);
diff --git a/OPF_Init.m b/OPF_Init.m
index 97142d6..f9eb03f 100644
--- a/OPF_Init.m
+++ b/OPF_Init.m
@@ -27,4 +27,5 @@ wD=ones(Busnum,1);
 %wD(Balance)=0;
 PD=1*PD0;
 %PD(PD==0)=.2;
+
 end
\ No newline at end of file