From 6285adc0d06fc20e5ee91936574b98cab7f54920 Mon Sep 17 00:00:00 2001 From: facat Date: Sat, 6 Apr 2013 21:14:59 +0800 Subject: [PATCH] =?UTF-8?q?=E5=88=A0=E9=99=A4=E4=BA=86=E4=B8=80=E4=BA=9B?= =?UTF-8?q?=E6=B2=A1=E7=94=A8=E7=9A=84=E6=96=87=E4=BB=B6=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- AssignXX.m | 23 ------ CalCost.m | 13 ---- DrawGap.asv | 6 -- DrawGap.m | 6 -- FormAA.asv | 12 --- FormAlphaD.m | 9 --- FormAlphaP.asv | 7 -- FormAlphaP.m | 9 --- FormG.asv | 26 ------- FormG.m | 34 -------- FormH.asv | 31 -------- FormH.m | 40 ---------- FormLw.asv | 16 ---- FormLw.m | 25 ------ FormLx.m | 4 - FormLxComa.m | 11 --- FormLz.asv | 22 ------ FormLz.m | 23 ------ FormYY.asv | 5 -- FormYY.m | 10 --- FormddgZW.m | 17 ---- Lineloss.asv | 28 ------- Modification.m | 28 ------- OPF.asv | 101 ------------------------ OPF.m | 121 ----------------------------- OPFWrongDataCheck.asv | 106 ------------------------- OPFWrongDataCheck.m | 106 ------------------------- OPF_Init.asv | 46 ----------- ObjectiveFun.asv | 19 ----- ReadMe.txt | 1 - Run_YALMIP.asv | 44 +++++++---- Run_YALMIP.m | 36 +++++---- SolveIt.asv | 42 ---------- SolveIt.m | 51 ------------ admmatrix.asv | 31 -------- func_ddf.m | 176 ------------------------------------------ func_ddg.asv | 14 ---- func_ddg.m | 23 ------ func_ddh.m | 58 -------------- func_deltF.asv | 19 ----- func_deltF.m | 15 ---- func_deltG.asv | 52 ------------- func_deltG.m | 58 -------------- func_deltH.asv | 8 -- func_deltH.m | 10 --- func_deltdeltF.asv | 11 --- func_deltdeltF.m | 11 --- jacobian_M.m | 20 ----- modifyadmmatrix.asv | 18 ----- modifyadmmatrix.m | 23 ------ pf.asv | 34 -------- test.m | 97 ----------------------- 52 files changed, 51 insertions(+), 1705 deletions(-) delete mode 100644 AssignXX.m delete mode 100644 CalCost.m delete mode 100644 DrawGap.asv delete mode 100644 DrawGap.m delete mode 100644 FormAA.asv delete mode 100644 FormAlphaD.m delete mode 100644 FormAlphaP.asv delete mode 100644 FormAlphaP.m delete mode 100644 FormG.asv delete mode 100644 FormG.m delete mode 100644 FormH.asv delete mode 100644 FormH.m delete mode 100644 FormLw.asv delete mode 100644 FormLw.m delete mode 100644 FormLx.m delete mode 100644 FormLxComa.m delete mode 100644 FormLz.asv delete mode 100644 FormLz.m delete mode 100644 FormYY.asv delete mode 100644 FormYY.m delete mode 100644 FormddgZW.m delete mode 100644 Lineloss.asv delete mode 100644 Modification.m delete mode 100644 OPF.asv delete mode 100644 OPF.m delete mode 100644 OPFWrongDataCheck.asv delete mode 100644 OPFWrongDataCheck.m delete mode 100644 OPF_Init.asv delete mode 100644 ObjectiveFun.asv delete mode 100644 ReadMe.txt delete mode 100644 SolveIt.asv delete mode 100644 SolveIt.m delete mode 100644 admmatrix.asv delete mode 100644 func_ddf.m delete mode 100644 func_ddg.asv delete mode 100644 func_ddg.m delete mode 100644 func_ddh.m delete mode 100644 func_deltF.asv delete mode 100644 func_deltF.m delete mode 100644 func_deltG.asv delete mode 100644 func_deltG.m delete mode 100644 func_deltH.asv delete mode 100644 func_deltH.m delete mode 100644 func_deltdeltF.asv delete mode 100644 func_deltdeltF.m delete mode 100644 jacobian_M.m delete mode 100644 modifyadmmatrix.asv delete mode 100644 modifyadmmatrix.m delete mode 100644 pf.asv delete mode 100644 test.m diff --git a/AssignXX.m b/AssignXX.m deleted file mode 100644 index 2ff0f36..0000000 --- a/AssignXX.m +++ /dev/null @@ -1,23 +0,0 @@ -function [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX1(XX,ContrlCount,RestraintCount,Busnum) -% deltX=XX(1:14); -% deltY=XX(15:24); -% deltZ=XX(25:38); -% deltW=XX(39:52); -% deltL=XX(53:66); -% deltU=XX(67:80); -deltX=XX(1:ContrlCount); -k1=ContrlCount+2*Busnum; -deltY=XX(ContrlCount+1:k1); -k2=k1; -k1=k2+RestraintCount; -deltZ=XX(k2+1:k1); -k2=k1; -k1=k2+RestraintCount; -deltW=XX(k2+1:k1); -k2=k1; -k1=k2+RestraintCount; -deltL=XX(k2+1:k1); -k2=k1; -k1=k2+RestraintCount; -deltU=XX(k2+1:k1); -end \ No newline at end of file diff --git a/CalCost.m b/CalCost.m deleted file mode 100644 index 1975144..0000000 --- a/CalCost.m +++ /dev/null @@ -1,13 +0,0 @@ -function CalCost(GenC,PG,PGi) -cost=GenC(:,2).*PG(PGi).^2+GenC(:,3).*PG(PGi)+GenC(:,4); -% Org_PG=[5; -% 2.5794]; -% book_PG=[5.5056; -% 2.1568]; -% cost2=GenC(:,2).*Org_PG(1:2).^2+GenC(:,3).*Org_PG(1:2)+GenC(:,4); -% cost3=GenC(:,2).*book_PG(1:2).^2+GenC(:,3).*book_PG(1:2)+GenC(:,4); -fprintf('总花费为%f\n',sum(cost,1)); -% fprintf('PF总花费为%f\n',sum(cost2,1)); -% fprintf('书上OPF总花费为%f\n',sum(cost3,1)); -% fprintf('较书上减少费用为为%f\n',sum(cost3,1)-sum(cost,1)); -end \ No newline at end of file diff --git a/DrawGap.asv b/DrawGap.asv deleted file mode 100644 index c9d0645..0000000 --- a/DrawGap.asv +++ /dev/null @@ -1,6 +0,0 @@ -function DrawGap(plotGap) -x=find(plotGap); -ts=size(x,2); - - -end \ No newline at end of file diff --git a/DrawGap.m b/DrawGap.m deleted file mode 100644 index f8d8727..0000000 --- a/DrawGap.m +++ /dev/null @@ -1,6 +0,0 @@ -function DrawGap(plotGap) -x=find(plotGap); -ts=size(x,2); -plot(1:ts,plotGap(1:ts)); - -end \ No newline at end of file diff --git a/FormAA.asv b/FormAA.asv deleted file mode 100644 index 0d82338..0000000 --- a/FormAA.asv +++ /dev/null @@ -1,12 +0,0 @@ -function AA=FormAA(L_1Z,deltG,U_1W,Hcoma,deltH) -tOnes=eye(14); -tZeros=zeros(14); -AA=[ - tOnes,L_1Z,tZeros,tZeros,tZeros,zeros(14,10); - tZeros,tOnes,tZeros,tZeros,-deltG',zeros(14,10); - tZeros,tZeros,tOnes,U_1W,tZeros,zeros(14,10); - tZeros,tZeros,tZeros,tOnes,deltG,zeros(14,10); - tZeros,tZeros,tZeros,tZeros,Hcoma,deltH; - zeros(10,14),zeros(10,14),zeros(10,14),zeros(10,14),deltH',zeros(10,10); - ]; -end \ No newline at end of file diff --git a/FormAlphaD.m b/FormAlphaD.m deleted file mode 100644 index 38816b9..0000000 --- a/FormAlphaD.m +++ /dev/null @@ -1,9 +0,0 @@ -function AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW) -tdeltZinx=find(deltZ<0); -tdeltWinx=find(deltW>0); -t1=-Init_Z(tdeltZinx)./deltZ(tdeltZinx)'; -t2=-Init_W(tdeltWinx)./deltW(tdeltWinx)'; -t3=[t1,t2]; -t4=min(t3); -AlphaD=0.9995*min([t4 1]); -end \ No newline at end of file diff --git a/FormAlphaP.asv b/FormAlphaP.asv deleted file mode 100644 index 3a8f5d9..0000000 --- a/FormAlphaP.asv +++ /dev/null @@ -1,7 +0,0 @@ -function AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU) -ti=deltL(delt) -t1=-Init_L./deltL'; -t2=-Init_U./deltU'; -t3=[t1,t2]; -t4=t3() -end \ No newline at end of file diff --git a/FormAlphaP.m b/FormAlphaP.m deleted file mode 100644 index a5757cc..0000000 --- a/FormAlphaP.m +++ /dev/null @@ -1,9 +0,0 @@ -function AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU) -tdeltLinx=find(deltL<0); -tdeltUinx=find(deltU<0); -t1=-Init_L(tdeltLinx)./deltL(tdeltLinx)'; -t2=-Init_U(tdeltUinx)./deltU(tdeltUinx)'; -t3=[t1,t2]; -t4=min(t3); -AlphaP=0.9995*min([t4 1]); -end \ No newline at end of file diff --git a/FormG.asv b/FormG.asv deleted file mode 100644 index eba25dc..0000000 --- a/FormG.asv +++ /dev/null @@ -1,26 +0,0 @@ -function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi) -PD_=PD(Loadi); -QD_=QD(Loadi); -Gsize=length(PGi)+length(PVi)+2*length(Loadi); -Mat_G=sdpvar(Gsize,1); -add=1; -for I=1:length(PGi) - Mat_G(add)=PG(I).value - add=add+1; -end -for I=1:length(PVi) - Mat_G(add)=QG(I).value - add=add+1; -end -for I=1:length(Loadi) - Mat_G(add)=QG(I).value - add=add+1; -end -% Mat_G=[ -% [PG(PGi).value]'; -% [QG(PVi).value]'; -% [PD(Loadi).value]'; -% [QD(Loadi).value]'; -% Volt; -% ]; -end \ No newline at end of file diff --git a/FormG.m b/FormG.m deleted file mode 100644 index 8cf6075..0000000 --- a/FormG.m +++ /dev/null @@ -1,34 +0,0 @@ -function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi) -PD_=PD(Loadi); -QD_=QD(Loadi); -Gsize=length(PGi)+length(PVi)+2*length(Loadi)+length(Volt); -Mat_G=sdpvar(Gsize,1); -add=1; -for I=1:length(PGi) - Mat_G(add)=PG(I).value; - add=add+1; -end -for I=1:length(PVi) - Mat_G(add)=QG(I).value; - add=add+1; -end -for I=1:length(Loadi) - Mat_G(add)=PD(I).value; - add=add+1; -end -for I=1:length(Loadi) - Mat_G(add)=QD(I).value; - add=add+1; -end -for I=1:length(Volt) - Mat_G(add)=Volt(I); - add=add+1; -end -% Mat_G=[ -% [PG(PGi).value]'; -% [QG(PVi).value]'; -% [PD(Loadi).value]'; -% [QD(Loadi).value]'; -% Volt; -% ]; -end \ No newline at end of file diff --git a/FormH.asv b/FormH.asv deleted file mode 100644 index d0226a8..0000000 --- a/FormH.asv +++ /dev/null @@ -1,31 +0,0 @@ -function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi,notLoadi,PGi,PVi,PD0,QD0) -%% -%QDcos=textread('D:\Project\最小化潮流\最小潮流算例\仙海919PDQDglys.txt'); -%QD(QD~=0)=PD(QD~=0)./tan(QDcos); -%QD(QD_NON_ZERO_IND)=QD_NON_ZERO; -%% -%PD(Loadi)=QD(Loadi)./tan(acos(0.98)); -AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); -Loadi=find(QD0~=0 | PD0~=0); -PDi=find(PD0~=0); -QDi=find(QD0~=0); -notLoadi=setdiff(1:Busnum,Loadi); -PGOnly=setdiff(PGi,PDi); -QGOnly=setdiff(PVi,QDi); -PDOnly=setdiff(PDi,PGi); -QDOnly=setdiff(QDi,PVi); -PGPD=intersect(PDi,PGi); -QGQD=intersect(QDi,PVi); -noPGPDQGQD=setdiff(1:Busnum,union(union(PDi,PGi),union(QDi,PVi))); -YCos=Y.*cos(AngleIJ); -YSin=Y.*sin(AngleIJ); -dP1=PG(ismember(PGi,PGOnly))-diag(Volt(PGOnly))*YCos(PGOnly,:)*Volt(PGOnly); -dQ1=QG(ismember(PVi,QGOnly))-diag(Volt(QGOnly))*YSin(QGOnly,:)*Volt(QGOnly); -dP2=PD(ismember(PDi,PDOnly))-diag(Volt(PDOnly))*YCos(PDOnly,:)*Volt(PDOnly); -dQ2=QD(ismember(QDi,QDOnly))-diag(Volt(QOnly))*YSin(QDOnly,:)*Volt(QDOnly); -dP3=PG(ismember(PGi,PGOnly))-diag(Volt(PGOnly))*YCos(PGOnly,:)*Volt(PGOnly); -dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt; - -Mat_H=[dP;dQ;]; - -end \ No newline at end of file diff --git a/FormH.m b/FormH.m deleted file mode 100644 index d885e01..0000000 --- a/FormH.m +++ /dev/null @@ -1,40 +0,0 @@ -function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi,notLoadi,PGi,PVi,PD0,QD0,dP,dQ) -%% -%QDcos=textread('D:\Project\最小化潮流\最小潮流算例\仙海919PDQDglys.txt'); -%QD(QD~=0)=PD(QD~=0)./tan(QDcos); -%QD(QD_NON_ZERO_IND)=QD_NON_ZERO; -%% -%PD(Loadi)=QD(Loadi)./tan(acos(0.98)); -AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); -Loadi=find(QD0~=0 | PD0~=0); -PDi=find(PD0~=0); -QDi=find(QD0~=0); -notLoadi=setdiff(1:Busnum,Loadi); -PGOnly=setdiff(PGi,PDi); -QGOnly=setdiff(PVi,QDi); -PDOnly=setdiff(PDi,PGi); -QDOnly=setdiff(QDi,PVi); -PGPD=intersect(PDi,PGi); -QGQD=intersect(QDi,PVi); -noPGPDQGQD=setdiff(1:Busnum,union(union(PDi,PGi),union(QDi,PVi))); -%% - YCos=Y.*cos(AngleIJ); - YSin=Y.*sin(AngleIJ); -% dP1=PG(ismember(PGi,PGOnly))-diag(Volt(PGOnly))*YCos(PGOnly,:)*Volt(PGOnly); -% dQ1=QG(ismember(PVi,QGOnly))-diag(Volt(QGOnly))*YSin(QGOnly,:)*Volt(QGOnly); -% dP2=PD(ismember(PDi,PDOnly))-diag(Volt(PDOnly))*YCos(PDOnly,:)*Volt(PDOnly); -% dQ2=QD(ismember(QDi,QDOnly))-diag(Volt(QOnly))*YSin(QDOnly,:)*Volt(QDOnly); -% dP3=PG(ismember(PGi,PGOnly))-diag(Volt(PGOnly))*YCos(PGOnly,:)*Volt(PGOnly); -% dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt; -%% - -for I=1:Busnum - dP(I)=PG(I).value-PD(I).value-Volt(I)*YCos(I,:)*Volt; - dQ(I)=QG(I).value-QD(I).value-Volt(I)*YSin(I,:)*Volt; -end -%% -% dP=[PG(:).value]'-[PD(:).value]'-diag(Volt)*YCos*Volt; -% dQ=[QG(:).value]'-[QD(:).value]'-diag(Volt)*YSin*Volt; -Mat_H=[dP;dQ;]; - -end \ No newline at end of file diff --git a/FormLw.asv b/FormLw.asv deleted file mode 100644 index 5dd6cf0..0000000 --- a/FormLw.asv +++ /dev/null @@ -1,16 +0,0 @@ -function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK) -%KK=999; -PU=1*GenU(:,2);%发电机有功上界 -QU=PVQU(:,1);%发电机无功上界 -%VoltU=(1.5+1/exp(KK))*ones(1,Busnum); -VoltU=1.5*ones(1,Busnum); -PDU=PD0(Loadi); -% PDU(PDU>0)=1200*PDU(PDU>0); -% PDU(PDU<0)=-800*PDU(PDU<0); -% PDU(PDU==0)=200; -PDU=1000000*ones(length(Loadi),1); -t1=([PU',QU',PDU',VoltU])'; -t2=Mat_G+Init_U'-t1; -Lw=t2; - -end \ No newline at end of file diff --git a/FormLw.m b/FormLw.m deleted file mode 100644 index ed01bcf..0000000 --- a/FormLw.m +++ /dev/null @@ -1,25 +0,0 @@ -function Lw=FormLw(Mat_G,GenU,Busnum,PVQU,PD0,QD0,Loadi,PF) -%PU=GenU(:,2);%发电机有功上界 -PU=8*ones(length(GenU(:,2)),1); -%QU=PVQU(:,1);%发电机无功上界 -QU=8*ones(length(PVQU(:,1)),1); -VoltU=(1.2)*ones(1,Busnum); -%VoltU=10*ones(1,Busnum); -PDU=PD0(Loadi); -PDU(PDU>0)=1.200*PDU(PDU>0); -PDU(PDU<0)=0.800*PDU(PDU<0); -PDU(PDU==0)=0.400; -%PDU=10*ones(length(Loadi),1); -QDU=QD0(Loadi); -QDU(QDU>0)=1.200*QDU(QDU>0); -QDU(QDU<0)=0.800*QDU(QDU<0); -QDU(QDU==0)=0.200; -%PF_=1.1*PF; -%PF_(PF_>1)=1; -%powerFactorU=2*log(1).*ones(length(Loadi),1); -%t1=([PU',QU',PDU',QDU',VoltU,powerFactorU'])'; -t1=([PU',QU',PDU',QDU',VoltU])'; -t2=Mat_G-t1; -Lw=t2; - -end \ No newline at end of file diff --git a/FormLx.m b/FormLx.m deleted file mode 100644 index 895f91e..0000000 --- a/FormLx.m +++ /dev/null @@ -1,4 +0,0 @@ -function Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W) -t1=deltF-deltH*Init_Y'-deltG*(Init_Z'+Init_W'); -Lx=t1; -end \ No newline at end of file diff --git a/FormLxComa.m b/FormLxComa.m deleted file mode 100644 index 9d25a12..0000000 --- a/FormLxComa.m +++ /dev/null @@ -1,11 +0,0 @@ -function LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx) -%t1=deltG*(Init_Z'+Init_W');%% -t2=Lul+diag(Init_Z)*Lz; -t3=inv(diag(Init_L)); -t4=t3*t2;% -t5=Luu-diag(Init_W)*Lw; -t6=inv(diag(Init_U)); -t7=t6*t5;% -t8=deltG*(t4+t7);%% -LxComa=Lx+t8; -end \ No newline at end of file diff --git a/FormLz.asv b/FormLz.asv deleted file mode 100644 index f8a1036..0000000 --- a/FormLz.asv +++ /dev/null @@ -1,22 +0,0 @@ -function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,QD0,Loadi,KK) -KK=999; -PL=GenL(:,2);%发电机有功下界 -QL=PVQL(:,1);%发电机无功下界 -VoltL=(0.8)*ones(1,Busnum); -%VoltL=-10*ones(1,Busnum); -PDL=PD0(Loadi); -PDL(PDL>0)=0.800*PDL(PDL>0); -PDL(PDL<0)=1.200*PDL(PDL<0); -PDL(PDL==0)=-0.400; -%PDL=-10*ones(length(Loadi),1); -% QDL=QD0(Loadi); -% QDL(QDL>0)=0.800*QDL(QDL>0); -% QDL(QDL<0)=1.200*QDL(QDL<0); -% QDL(QDL==0)=-0.200; -powerFactorL=2*log(cos(0.90))*ones(length(Loadi),1); -%t1=([PL',QL',PDL',QDL',VoltL,powerFactorL'])'; -t1=([PL',QL',PDL',QDL,VoltL,powerFactorL'])'; -t2=Mat_G-Init_L'-t1; -Lz=t2; - -end \ No newline at end of file diff --git a/FormLz.m b/FormLz.m deleted file mode 100644 index e7ea0a7..0000000 --- a/FormLz.m +++ /dev/null @@ -1,23 +0,0 @@ -function Lz=FormLz(Mat_G,GenL,Busnum,PVQL,PD0,QD0,Loadi,PF) -%PL=GenL(:,2);%发电机有功下界 -PL=-60*ones(length(GenL(:,2)),1); -%QL=PVQL(:,1);%发电机无功下界 -QL=-60*ones(length(PVQL(:,1)),1); -VoltL=(0.8)*ones(1,Busnum); -%VoltL=-10*ones(1,Busnum); -PDL=PD0(Loadi); -PDL(PDL>0)=0.800*PDL(PDL>0); -PDL(PDL<0)=1.200*PDL(PDL<0); -PDL(PDL==0)=-0.400; -%PDL=-10*ones(length(Loadi),1); -QDL=QD0(Loadi); -QDL(QDL>0)=0.800*QDL(QDL>0); -QDL(QDL<0)=1.200*QDL(QDL<0); -QDL(QDL==0)=-0.200; -%powerFactorL=2*log(0.75).*ones(length(Loadi),1); -%t1=([PL',QL',PDL',QDL',VoltL,powerFactorL'])'; -t1=([PL',QL',PDL',QDL',VoltL])'; -t2=Mat_G-t1; -Lz=t2; - -end \ No newline at end of file diff --git a/FormYY.asv b/FormYY.asv deleted file mode 100644 index 85d2a13..0000000 --- a/FormYY.asv +++ /dev/null @@ -1,5 +0,0 @@ -function FormYY(Init_L,Lul.Lz,Init_U,Luu,Lw,Lz,LxComa) -t=[ - - -] -end \ No newline at end of file diff --git a/FormYY.m b/FormYY.m deleted file mode 100644 index cbad106..0000000 --- a/FormYY.m +++ /dev/null @@ -1,10 +0,0 @@ -function YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx) -YY=[ - Lx; - -Ly; - -Lz; - -Lw; - -Lul; - -Luu; -]; -end \ No newline at end of file diff --git a/FormddgZW.m b/FormddgZW.m deleted file mode 100644 index c68abf9..0000000 --- a/FormddgZW.m +++ /dev/null @@ -1,17 +0,0 @@ -function [ ddgZW ] = FormddgZW(PGi,PVi,Busnum,ContrlCount,Loadi,PD,QD,Init_Z,Init_W) - -ddg1=sparse(size(PVi,1)+size(PGi,1),ContrlCount); -%计算W和Z应该乘的位置 -rePos1=length(PGi)+length(PVi)+length(Loadi)*2+Busnum+1; -rePos2=rePos1+length(Loadi)-1; -t1=sparse(length(Loadi),size(PVi,1)+size(PGi,1)); -t2=diag(-2./(PD(Loadi).^2)-2*(QD(Loadi).^2-PD(Loadi).^2)./(PD(Loadi).^2+QD(Loadi).^2))*diag(Init_Z(rePos1:rePos2)+Init_W(rePos1:rePos2)); -t3=4*diag(PD(Loadi).*QD(Loadi)./((PD(Loadi).^2+QD(Loadi).^2).^2))*diag(Init_Z(rePos1:rePos2)+Init_W(rePos1:rePos2)); -t4=sparse(length(Loadi),2*Busnum); -ddg2=[t1,t2,t3,t4]; -t2=diag(-2*(PD(Loadi).^2-QD(Loadi).^2)./(PD(Loadi).^2+QD(Loadi).^2))*diag(Init_Z(rePos1:rePos2)+Init_W(rePos1:rePos2)); -ddg3=[t1,t3,t2,t4]; -ddg4=sparse(Busnum*2,ContrlCount); -ddgZW=[ddg1;ddg2;ddg3;ddg4]; -end - diff --git a/Lineloss.asv b/Lineloss.asv deleted file mode 100644 index 3132a3d..0000000 --- a/Lineloss.asv +++ /dev/null @@ -1,28 +0,0 @@ -%% 计算线损 -function Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,k0,Volt,Angle) -%format long -fprintf('功率为有名值\n'); -fprintf('节点号\t节点号\t有功') -cmpVolt=Volt.*cos(Angle)+1i*Volt.*sin(Angle); -cmpVolt=cmpVolt'; -y0=1i*Lineb2; -%yj0=1i*standardinput(:,7); -yij=1./(Liner+1i*Linex); -%% 线路损耗 -Sij=cmpVolt(Linei).*conj( cmpVolt(Linei) .* y0 + ( cmpVolt(Linei)- cmpVolt(Linej) ).*yij ); -Sji=cmpVolt(Linej).*conj( cmpVolt(Linej) .*y0 + ( cmpVolt(Linej)- cmpVolt(Linei) ).*yij ); -%Sij1==Sij2 -deltLineS=Sij+Sji; -dispLineloss=[Linei Linej real(deltLineS)*100 imag(deltLineS)*100]; -%full(dispLineloss) -dispLineloss=sortrows(dispLineloss,-3); -full(dispLineloss) -%% 以下是变压器损耗 -yij=1./(Transr+1i*Transx); -Sij=cmpVolt(Transi)./k0.*conj( ( cmpVolt(Transi)./k0- cmpVolt(Transj) ).*yij ); -Sji=cmpVolt(Transj).*conj( ( cmpVolt(Transj)- cmpVolt(Transi)./k0 ).*yij ); -deltTransS=Sij+Sji; -dispTransloss=[Transi Transj real(deltTransS)*100 imag(deltTransS)*100]; -dispTransloss=sortrows(dispTransloss,-3); -full(dispTransloss) -end \ No newline at end of file diff --git a/Modification.m b/Modification.m deleted file mode 100644 index 7da8912..0000000 --- a/Modification.m +++ /dev/null @@ -1,28 +0,0 @@ -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); -%fprintf('AlphaD %f\n',full(AlphaD)); - -Init_Z=Init_Z+AlphaD*deltZ'; -Init_L=Init_L+AlphaP*deltL'; -Init_W=Init_W+AlphaD*deltW'; -Init_U=Init_U+AlphaP*deltU'; -Init_Y=Init_Y+AlphaD*deltY'; -%PG(PGi)=PG(PGi)+deltX(1:size(PGi,1)); -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)*2); -PD(Loadi)=PD(Loadi)+AlphaP*t(1:length(Loadi)); -%PD(PD<0)=-PD(PD<0); -QD(Loadi)=QD(Loadi)+AlphaP*t(length(Loadi)+1:length(Loadi)*2); -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 -balVolt=Volt(Balance); -Volt=Volt+AlphaP*t(1:Busnum); -Volt(Balance)=balVolt; -UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum); -end \ No newline at end of file diff --git a/OPF.asv b/OPF.asv deleted file mode 100644 index 10abed9..0000000 --- a/OPF.asv +++ /dev/null @@ -1,101 +0,0 @@ -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('ieee4.dat'); -%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); -%pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); -%pf('c:/file31.txt'); - - -%% 计算功率因数 -%atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0)); -Volt; -UAngel*180/3.1415926; -%% 通过潮流计算PG -AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); -PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt'; -QGBal=diag(Volt)*Y.*sin(AngleIJ)*Volt'; -%% 初值-即测量值 -PG0=PG; -QG0=QG; -PD0=PD; -QD0=QD; -PDReal=PD;%真值 -%PD0(12)=PD0(12)+0.001; -%% -PG0(Balance)=PGBal(Balance); -QG0(Balance)=QGBal(Balance); -%% -[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD); -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)*2+Busnum*2; -kmax=60; -%% 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; - end - plotGap(KK+1)=Gap; - Init_u=Gap/2/RestraintCount*CenterA; - AngleIJMat=0; - %% 开始计算OPF - %% 形成等式约束的雅克比 - deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi); - %% 形成不等式约束的雅克比 - deltG=func_deltG(Busnum,PVi,PGi,Loadi); - %% - L_1Z=diag(Init_Z./Init_L); - U_1W=diag(Init_W./Init_U); - %% 形成海森阵 - deltdeltF=func_deltdeltF(PVi,wPG,wQG,wD,ContrlCount); - %% 形成ddHy - ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount); - %% 开始构建ddg - ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi); - %% 开始构建deltF - deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wD,PG0,QG0,PD0,PD,Busnum,Loadi); - - %% 形成方程矩阵 - 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,QD,Loadi); - Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND); - Ly=Mat_H; - Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,QD0,Loadi,KK); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,QD0,Loadi,KK); - Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); - YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); - %% 开始解方程 - fprintf(') - 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,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 -fprintf('迭代次数%d\n',KK); -%fprintf('目标值%f\n',full(ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi))); -DrawGap(plotGap); -%% -%Volt=full(Volt'); -%PD=full(PD); -%% 统计PD误差 -% absPDLoad=abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) ); -absPDLoad=abs( (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi) ); -maxPDError=max(absPDLoad); -disp('index'); -Loadi(absPDLoad==maxPDError); -%% 计算总线损 -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.m b/OPF.m deleted file mode 100644 index 9c52595..0000000 --- a/OPF.m +++ /dev/null @@ -1,121 +0,0 @@ -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('ieee14.dat'); -%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); -%pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); -%pf('c:/file31.txt'); - -Loadi=QD~=0 | PD~=0; -%% 计算功率因数 -PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2)) -Volt; -UAngel*180/3.1415926; -%% 通过潮流计算PG -AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); -PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt'; -QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt'; -%% 初值-即测量值 -PG0=PG; -QG0=QG; -PD0=PD; -QD0=QD; -PDReal=PD;%真值 -QDReal=QD;%真值 -%PD0(12)=PD0(12)+0.001; -%% -PG0(Balance)=PGBal(Balance); -QG0(Balance)=QGBal(Balance); -QG0(PVi)=QGBal(PVi); -PG(Balance)=PGBal(Balance); -QG(PVi)=QGBal(PVi); -%% -AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); -dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt'; -dPD=abs(dP./PD); -dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt'; -dQD=abs(dQ./QD); -maxdPQ=max([dPD(dPD<10);dQD(dQD<10)]) - -%% -[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD); -Gap=(Init_L*Init_Z'-Init_U*Init_W'); -KK=0; -plotGap=zeros(1,60); -ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+Busnum*2; -kmax=150; -%% 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; - end - plotGap(KK+1)=Gap; - Init_u=Gap/2/RestraintCount*CenterA; - AngleIJMat=0; - %% 开始计算OPF - %% 形成等式约束的雅克比 - deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi); - %% 形成不等式约束的雅克比 - deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD); - %% - L_1Z=diag(Init_Z./Init_L); - U_1W=diag(Init_W./Init_U); - %% 形成海森阵 - deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount); - %% 形成ddHy - ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount); - %% 开始构建ddg - ddg=func_ddg(PGi,PVi,Busnum,ContrlCount,RestraintCount,Loadi,PD,QD,Init_Z,Init_W); - %% 开始构建deltF - deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD0,QD,Busnum,Loadi); - - %% 形成方程矩阵 - 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,QD,Loadi); - Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND,Loadi); - Ly=Mat_H; - Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,QD0,Loadi,KK,PF); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,QD0,Loadi,KK,PF); - Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); - YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); - %% 开始解方程 - fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1)); - 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,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 -fprintf('迭代次数%d\n',KK); -fprintf('目标值%f\n',full(ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi))); -DrawGap(plotGap); -%% -%Volt=full(Volt'); -%PD=full(PD); -%% 统计PD误差 -% absPDLoad=abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) ); -absPDLoad=abs( (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi) ); -absQDLoad=abs( (QD(Loadi)-QDReal(Loadi))./QDReal(Loadi) ); -maxPDError=max([absPDLoad]) -maxQDError=max([absQDLoad]) -disp('index'); -Loadi(absPDLoad==maxPDError); -%% 计算总线损 -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 -AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); -dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt'; -dPD=abs(dP./PD); -dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt'; -dQD=abs(dQ./QD); -maxdPQ=max([dPD(dPD<10);dQD(dQD<10)]) -PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2)) \ No newline at end of file diff --git a/OPFWrongDataCheck.asv b/OPFWrongDataCheck.asv deleted file mode 100644 index 04bce85..0000000 --- a/OPFWrongDataCheck.asv +++ /dev/null @@ -1,106 +0,0 @@ -clc - - clear - for II=1:16 - tic - [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('D:\Project\最小化潮流\最小潮流算例\金湖924（2-1）_0.5_85%.txt'); - %pf('c:/file31.txt'); - %pf('ieee118PG.dat'); - - %% 计算功率因数 - atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0)); - Volt; - UAngel*180/3.1415926; - %% 通过潮流计算PG - AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); - PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt'; - - %% 初值-即测量值 - PG0=PG; - PD0=PD; - - PDReal=PD;%真值 - - %% - PG0(Balance)=PGBal(Balance); - %% - [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,PD); - PD0(Loadi(II))=PD0(Loadi(II))*(1+0.018); - 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; - kmax=60; - %% 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; - end - plotGap(KK+1)=Gap; - Init_u=Gap/2/RestraintCount*CenterA; - AngleIJMat=0; - %% 开始计算OPF - %% 形成等式约束的雅克比 - deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi); - %% 形成不等式约束的雅克比 - deltG=func_deltG(Busnum,PVi,PGi,Loadi); - %% - L_1Z=diag(Init_Z./Init_L); - U_1W=diag(Init_W./Init_U); - %% 形成海森阵 - deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount); - %% 形成ddHy - ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount); - %% 开始构建ddg - ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi); - %% 开始构建deltF - deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi); - - %% 形成方程矩阵 - 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,Loadi); - Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND); - Ly=Mat_H; - Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK); - Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); - YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); - %% 开始解方程 - 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); - Gap=(Init_L*Init_Z'-Init_U*Init_W'); - KK=KK+1; - end - fprintf('迭代次数%d\n',KK); - ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi); - %DrawGap(plotGap); - %% - %Volt=full(Volt'); - %PD=full(PD); - %% 统计PD误差 - % absPDLoad=abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) ); - absPDLoad=abs( (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi) ); - maxPDError=max(absPDLoad); - %disp('index') - LoadiArray=Loadi(absPDLoad==maxPDError); - if length(LoadiArray)>1 - disp('没找出') - Loadi(II) - end - if length(LoadiArray)==1 - if LoadiArray~=Loadi(II) - disp('没找出') - Loadi(II) - end - end - toc; - -end - diff --git a/OPFWrongDataCheck.m b/OPFWrongDataCheck.m deleted file mode 100644 index 16b6b1a..0000000 --- a/OPFWrongDataCheck.m +++ /dev/null @@ -1,106 +0,0 @@ -clc - - clear - for II=1:53 - tic - [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('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); - %pf('c:/file31.txt'); - %pf('ieee118PG.dat'); - - %% 计算功率因数 - atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0)); - Volt; - UAngel*180/3.1415926; - %% 通过潮流计算PG - AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); - PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt'; - - %% 初值-即测量值 - PG0=PG; - PD0=PD; - - PDReal=PD;%真值 - - %% - PG0(Balance)=PGBal(Balance); - %% - [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,PD); - PD0(Loadi(II))=PD0(Loadi(II))*(1+0.086); - 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; - kmax=60; - %% 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; - end - plotGap(KK+1)=Gap; - Init_u=Gap/2/RestraintCount*CenterA; - AngleIJMat=0; - %% 开始计算OPF - %% 形成等式约束的雅克比 - deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi); - %% 形成不等式约束的雅克比 - deltG=func_deltG(Busnum,PVi,PGi,Loadi); - %% - L_1Z=diag(Init_Z./Init_L); - U_1W=diag(Init_W./Init_U); - %% 形成海森阵 - deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount); - %% 形成ddHy - ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount); - %% 开始构建ddg - ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi); - %% 开始构建deltF - deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi); - - %% 形成方程矩阵 - 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,Loadi); - Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND); - Ly=Mat_H; - Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK); - Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); - YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); - %% 开始解方程 - 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); - Gap=(Init_L*Init_Z'-Init_U*Init_W'); - KK=KK+1; - end - fprintf('迭代次数%d\n',KK); - ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi); - %DrawGap(plotGap); - %% - %Volt=full(Volt'); - %PD=full(PD); - %% 统计PD误差 - % absPDLoad=abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) ); - absPDLoad=abs( (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi) ); - maxPDError=max(absPDLoad); - %disp('index') - LoadiArray=Loadi(absPDLoad==maxPDError); - if length(LoadiArray)>1 - disp('没找出') - Loadi(II) - end - if length(LoadiArray)==1 - if LoadiArray~=Loadi(II) - disp('没找出') - Loadi(II) - end - end - toc; - -end - diff --git a/OPF_Init.asv b/OPF_Init.asv deleted file mode 100644 index 2eb039a..0000000 --- a/OPF_Init.asv +++ /dev/null @@ -1,46 +0,0 @@ -function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD) -Loadi=find(QD~=0 | PD~=0); -PDi=find(PD~=0); -QDi=find(QD~=0); -notLoadi=setxor(1:Busnum,Loadi); -PGOnly=setdiff(PGi,PDi); -QGOnly=setdiff(PVi,QDi); -PDOnly=setdiff(PDi,PGi); -QDOnly=setdiff(QDi,PVi); -PGPD=intersec(PDi,PGi); -QGQD=intersec(QDi,PVi); -noPGPDQGQD=setdiff(1:B) -%Loadi=[1:Busnum]'; -RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*1+Busnum*1; %约束条件数,放开所有QD,不单独约束QD上下限 -RestraintCount=RestraintCount+length(Loadi); %加上功率因数，用功率因数约束QD -t_Bal_volt=Volt(Balance); -Volt=sparse(1*ones(1,Busnum)); -Volt(Balance)=t_Bal_volt; -UAngel=sparse(1,Busnum); -Init_Z=sparse(ones(1,RestraintCount)); -Init_W=sparse(-1*ones(1,RestraintCount)); -Init_L=1*sparse(ones(1,RestraintCount)); -Init_U=1*sparse(ones(1,RestraintCount)); -Init_Y=sparse(1,2*Busnum);%与学姐一致 -tPU=sparse(GenU(:,2));% 发电机有功上限 -tQU=sparse(PVQU(:,1));% 无功上限 -tPL=sparse(GenL(:,2));% 发电机有功下限 -tQL=sparse(PVQL(:,1));% 无功下限 -PG(PGi)=(tPU+tPL)/2; -QG(PVi)=(tQU+tQL)/2; -wPG=1*ones(size(PGi,1),1); -wQG=1*ones(size(PVi,1),1); -%randInt=randperm(size(Loadi,1)); -%randPDind=randInt(1:10); -randPDind=0; -wPD=1*ones(size(Loadi,1),1); -wQD=1.5*ones(size(Loadi,1),1); -%wD(randPDind)=0;%一些负荷不约束 -%wD(7)=0; -% wD(11)=0; -PD=1*PD0; -%powerFacter=0.98; -%QD=PD*sqrt((1-powerFacter^2)/powerFacter^2); -QD=QD0; - -end \ No newline at end of file diff --git a/ObjectiveFun.asv b/ObjectiveFun.asv deleted file mode 100644 index a892fa0..0000000 --- a/ObjectiveFun.asv +++ /dev/null @@ -1,19 +0,0 @@ -function [out_arg]=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi) -% t1=[PG(PGi).value]'-PG0(PGi); -% t2=wPG.*t1.*t1; -t2=sdpvar(1,1); -for I=1:length(PGi) - t2=t2+wPG(I)*(PG(I).value-PG0(I)); -end - -% t1=[QG(PVi).value]'-QG0(PVi); -% t3=wQG.*t1.*t1; -t3=sdpvar(1,1); -for I=1:length(PVi) - t3=t3+wPG(I)*(QG(I).value-QG0(I)); -end -% t4=wPD.*(([PD(Loadi).value]'-PD0(Loadi)).*([PD(Loadi).value]'-PD0(Loadi))); -for I=1:length(L) -% t5=wQD.*(([QD(Loadi).value]'-QD0(Loadi)).*([QD(Loadi).value]'-QD0(Loadi))); -out_arg= sum(t2)+sum(t3)+sum(t4)+sum(t5); -end \ No newline at end of file diff --git a/ReadMe.txt b/ReadMe.txt deleted file mode 100644 index 22b1747..0000000 --- a/ReadMe.txt +++ /dev/null @@ -1 +0,0 @@ -对照学姐给的公式 \ No newline at end of file diff --git a/Run_YALMIP.asv b/Run_YALMIP.asv index b691b8d..628fdb5 100644 --- a/Run_YALMIP.asv +++ b/Run_YALMIP.asv @@ -5,7 +5,7 @@ tic [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('ieee4.dat'); + pf('ieee30.dat'); %% 潮流等式 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); @@ -37,40 +37,52 @@ maxdPQ=max([dPD(dPD<10);dQD(dQD<10)]); %% 定义变量 Volt=sdpvar(Busnum,1); UAngel=sdpvar(Busnum,1); -PG=sdpvar(Busnum,1); -QG=sdpvar(Busnum,1); +% PG=sdpvar(Busnum,1); +% QG=sdpvar(Busnum,1); PD=sdpvar(Busnum,1); QD=sdpvar(Busnum,1); -%AngleIJ=sdpvar(Busnum,Busnum,'full'); +AngleIJ=sdpvar(Busnum,Busnum,'full'); %% 目标函数 Objective=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi); -AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); +%AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); +%% 赋初值，可以加快求解速度。 +assign(Volt(:),1); +assign(UAngel(:),0); +assign(PD(:),PD0(:)); +assign(QD(:),QD0(:)); +% assign(PG(:),PG0(:)); +% assign(QG(:),QG0(:)); %% YALMIP部分 -Constraints = [%AngleIJ==sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum), ... - AngleIJ==0 - PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt==0, ... - QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt==0, ... +dP=PG0-PD-diag(Volt)*Y.*cos( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; +dQ=QG0-QD-diag(Volt)*Y.*sin( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; +Loadi=PD0~=0 | QD0~=0; +Constraints = [%AngleIJ-sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum)==0, ... + dP==0, ... + dQ==0, ... PG(setxor(1:Busnum, PGi) )==0, ... QG(setxor(1:Busnum, PVi) )==0, ... - PD(setdiff(1:Busnum,PD0~=0))==0, ... - QD(setdiff(1:Busnum,QD0~=0))==0, ... + PD(PD0==0)==0, ... + QD(QD0==0)==0, ... 0.9*ones(Busnum,1)<=Volt<=1.1*ones(Busnum,1), ... - GenL(:,2)<=PG(PGi)<=GenU(:,2), ... - PVQL<=QG(PVi)<=PVQU, ... + -60*ones(Busnum,1)<=PG<=60*ones(Busnum,1), ... + -60*ones(Busnum,1)<=QG<=60*ones(Busnum,1) ]; -options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt'); +options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt','usex0','1'); sol = solvesdp(Constraints,Objective,options); if sol.problem == 0 fprintf('Volt\n'); - dvolt=double(Volt); + dvolt=double(Volt) fprintf('VoltAngle\n'); - dVangle=double(UAngel); + dVangle=double(UAngel) fprintf('ojb\n'); optimalObj=double(Objective) + double(PG)-PG0; sol else display('Hmm, something went wrong!'); sol.info + sol.solveroutput yalmiperror(sol.problem) end + toc diff --git a/Run_YALMIP.m b/Run_YALMIP.m index be30f0e..628fdb5 100644 --- a/Run_YALMIP.m +++ b/Run_YALMIP.m @@ -37,42 +37,52 @@ maxdPQ=max([dPD(dPD<10);dQD(dQD<10)]); %% 定义变量 Volt=sdpvar(Busnum,1); UAngel=sdpvar(Busnum,1); -PG=sdpvar(Busnum,1); -QG=sdpvar(Busnum,1); +% PG=sdpvar(Busnum,1); +% QG=sdpvar(Busnum,1); PD=sdpvar(Busnum,1); QD=sdpvar(Busnum,1); AngleIJ=sdpvar(Busnum,Busnum,'full'); %% 目标函数 Objective=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi); %AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); +%% 赋初值，可以加快求解速度。 +assign(Volt(:),1); +assign(UAngel(:),0); +assign(PD(:),PD0(:)); +assign(QD(:),QD0(:)); +% assign(PG(:),PG0(:)); +% assign(QG(:),QG0(:)); %% YALMIP部分 -Constraints = [AngleIJ==sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum), ... - PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt==0, ... - QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt==0, ... +dP=PG0-PD-diag(Volt)*Y.*cos( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; +dQ=QG0-QD-diag(Volt)*Y.*sin( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; +Loadi=PD0~=0 | QD0~=0; +Constraints = [%AngleIJ-sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum)==0, ... + dP==0, ... + dQ==0, ... PG(setxor(1:Busnum, PGi) )==0, ... QG(setxor(1:Busnum, PVi) )==0, ... PD(PD0==0)==0, ... QD(QD0==0)==0, ... 0.9*ones(Busnum,1)<=Volt<=1.1*ones(Busnum,1), ... - %GenL(:,2)<=PG(PGi)<=GenU(:,2), ... - %PVQL<=QG(PVi)<=PVQU, ... - -20*ones(Busnum,1)<=PG<=20*ones(Busnum,1), ... - -20*ones(Busnum,1)<=QG<=20*ones(Busnum,1) + -60*ones(Busnum,1)<=PG<=60*ones(Busnum,1), ... + -60*ones(Busnum,1)<=QG<=60*ones(Busnum,1) ]; -options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt'); +options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt','usex0','1'); sol = solvesdp(Constraints,Objective,options); if sol.problem == 0 fprintf('Volt\n'); - dvolt=double(Volt); + dvolt=double(Volt) fprintf('VoltAngle\n'); - dVangle=double(UAngel); + dVangle=double(UAngel) fprintf('ojb\n'); optimalObj=double(Objective) - double(PG)-PG0 + double(PG)-PG0; sol else display('Hmm, something went wrong!'); sol.info + sol.solveroutput yalmiperror(sol.problem) end + toc diff --git a/SolveIt.asv b/SolveIt.asv deleted file mode 100644 index dceba29..0000000 --- a/SolveIt.asv +++ /dev/null @@ -1,42 +0,0 @@ -function 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) -LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx); -H=-deltdeltF+ddh;%+ddg*(Init_Z'+Init_W'); -t1=diag(Init_L.\Init_Z-Init_U.\Init_W); -t2=-deltG*( t1 )*deltG'; -aa=[ - (H+t2),deltH; - deltH',zeros(size(Init_Y,2)); - ]; -yy=[LxComa;-Ly]; -%% 平衡节点电压不变 -t=size(PVi,1)+size(PGi,1)+size(Loadi,1); -aa(t+Balance,:)=0; -aa(:,t+Balance)=0; -%aa(t+Balance,t+Balance)=1; -aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance)),ContrlCount+2*Busnum,ContrlCount+2*Busnum); -deltG(t+Balance,:)=0; -%% -t=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1; -aa(t+Balance,:)=0; -aa(:,t+Balance)=0; -%aa(t+Balance,t+Balance)=1; -aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance)),ContrlCount,ContrlCount); -deltG(t+Balance,:)=0; -%% -dxdy=aa\yy; -dX=dxdy(1:ContrlCount); -dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum); -dL=Lz+deltG'*dX; -dU=-Lw-deltG'*dX; -dZ=-diag(Init_L)\Lul-diag(Init_L)\diag(Init_Z)*dL; -dW=-diag(Init_U)\Luu-diag(Init_U)\diag(Init_W)*dU; -XX=[ - dX; - dY; - dZ; - dW; - dL; - dU; - -]; -end \ No newline at end of file diff --git a/SolveIt.m b/SolveIt.m deleted file mode 100644 index a8c7fdc..0000000 --- a/SolveIt.m +++ /dev/null @@ -1,51 +0,0 @@ -function 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) -LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx); -H=-deltdeltF+ddh+ddg; -%t1=diag(Init_L.\Init_Z-Init_U.\Init_W); -t1=diag(Init_Z./Init_L-Init_W./Init_U); -t2=-deltG*( t1 )*deltG'; -aa=[ - (H+t2),deltH; - deltH',zeros(size(Init_Y,2)); - ]; -yy=[LxComa;-Ly]; - -%% 平衡节点电压不变 -t=size(PVi,1)+size(PGi,1)+size(Loadi,1)*2; -yy(t+Balance)=0; -aa(t+Balance,:)=0; -aa(:,t+Balance)=0; -%aa(t+Balance,t+Balance)=1; -aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); -deltG(t+Balance,:)=0; -%% -t=size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+Busnum*1; -yy(t+Balance)=0; -aa(t+Balance,:)=0; -aa(:,t+Balance)=0; -%aa(t+Balance,t+Balance)=1; -aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); -deltG(t+Balance,:)=0; -%% - -dxdy=aa\yy; -%% KLU -%spy(aa) -%dxdy = klu(aa,'\',full(yy)); -%% -dX=dxdy(1:ContrlCount); -dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum); -dL=Lz+deltG'*dX; -dU=-Lw-deltG'*dX; -dZ=-diag(Init_L)\Lul-diag(Init_L)\diag(Init_Z)*dL; -dW=-diag(Init_U)\Luu-diag(Init_U)\diag(Init_W)*dU; -XX=[ - dX; - dY; - dZ; - dW; - dL; - dU; - -]; -end \ No newline at end of file diff --git a/admmatrix.asv b/admmatrix.asv deleted file mode 100644 index 1896345..0000000 --- a/admmatrix.asv +++ /dev/null @@ -1,31 +0,0 @@ -function [G,B,GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori... - ,Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb) -%************************************************************************** -% 程序功能 : 子函数——形成节点导纳矩阵Y -% 编 者： -% 编制时间：2010.12 -%************************************************************************** -%% 支路导纳计算 -G = -sparse(Linei,Linej,Liner./(Liner.^2+Linex.^2),Busnum,Busnum) - sparse(Linej,Linei,Liner./(Liner.^2+Linex.^2),Busnum,Busnum); -G = G - sparse(1:Busnum,1:Busnum,sum(G,2)'); % 计算各线路支路电导 -B = sparse(Linei,Linej,Linex./(Liner.^2+Linex.^2),Busnum,Busnum)+sparse(Linej,Linei,Linex./(Liner.^2+Linex.^2),Busnum,Busnum); -B = B - sparse(1:Busnum,1:Busnum,sum(B,2)')+sparse(Linei,Linei,Lineb,Busnum,Busnum)+sparse(Linej,Linej,Lineb,Busnum,Busnum); -%% 变压器支路计算 -if Transfori>0 - mr = Transforr./(Transforr.^2+Transforx.^2); % 计算变压器支路电导 - mx = -Transforx./(Transforr.^2+Transforx.^2); % 计算变压器支路电纳 - G = G-sparse(Transfori,Transforj,mr./Transfork0,Busnum,Busnum)-sparse(Transforj,Transfori,mr./Transfork0,Busnum,Busnum)... - +sparse(Transfori,Transfori,mr./Transfork0./Transfork0,Busnum,Busnum)+sparse(Transforj,Transforj,mr,Busnum,Busnum); - B = B-sparse(Transfori,Transforj,mx./Transfork0,Busnum,Busnum)-sparse(Transforj,Transfori,mx./Transfork0,Busnum,Busnum)... - +sparse(Transfori,Transfori,mx./Transfork0./Transfork0,Busnum,Busnum)+sparse(Transforj,Transforj,mx,Busnum,Busnum); -end -%% 接地支路计算 -if Branchi>0 % 判断有无接地支路 - B = B+sparse(Branchi,Branchi,Branchb,Busnum,Busnum); -end -%% 化作极坐标形式 -GB = G+B.*1i; %将电导，电纳合并，写成复数形式 -Y = abs(GB); %求节点导纳幅值 -[r,c] = find(Y); -Angle = angle(GB(GB~=0)); %求节点导纳角度 -%Angle=angle(GB); \ No newline at end of file diff --git a/func_ddf.m b/func_ddf.m deleted file mode 100644 index 594161d..0000000 --- a/func_ddf.m +++ /dev/null @@ -1,176 +0,0 @@ -function ddf=func_ddh(AngleIJMat,GB,Volt,Init_Y,Busnum) -%% deltaPi/deltaThytai_deltaThytaj 非对角元素 -t1=-Volt'*Volt; -t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -dPdTidTj=t1.*t2; %%(保留了对角元素的) -dPidTidTj_2=dPdTidTj-diag(diag(dPdTidTj));%去掉了对角元素的 -t3=repmat(Init_Y',1,size(Init_Y,2));%乘y的系数 -dPidTidTj_2=dPidTidTj_2.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2)); -t2=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat); -t3=t1.*t2; -dPjdTidTj=t3-diag(diag(t3));%去掉了对角元素的 -t3=repmat(Init_Y,size(Init_Y,2),1); -dPjdTidTj=dPjdTidTj.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数 -dPdTidTj=dPidTidTj_2+dPjdTidTj;%最终非对角元素 @@ -%% deltaP/deltaThyta_deltaThyta 对角元素 -t1=sum(-dPidTidTj_2,2); -t2=diag(t1'.*Init_Y(1:2:size(Init_Y,2)));%乘y的系数 -t3=sum(-dPidTidTj_2,1); -t4=diag(t3.*Init_Y(1:2:size(Init_Y,2)));%乘y的系数 -dPdTidTi=t2+t4;%%最终对角元素 @@ -%% deltaP/deltaThytai_dVi 对角元素 -t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -t2=diag(t1); -t3=t1-diag(t2);%去掉了对角元素的 -t4=sum(t3,2); -t4=t4'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数 -dPidTidVi=diag(t4); -dPidTjdVj=-t3; -t6=sum(dPidTjdVj,1);%乘y的系数 -t6=t6.*Init_Y(1:2:size(Init_Y,2));%乘y的系数 -dPdTidVi=dPidTidVi+diag(t6);%%最终对角元素 @@ -%% deltaP/deltaThytai_dVj 非对角元素 -t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -dPidTidVj=t1-diag(diag(t1));%%去掉对角元素的 -t2=repmat(Init_Y',1,size(Init_Y,2)); -dPidTidVj=dPidTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数 -t1=ones(Busnum,1)*Volt.*real(GB).*(sin(AngleIJMat)+imag(GB).*cos(AngleIJMat)); -dPjdTidVj=t1-diag(diag(t1));%去掉对角元素的 -t2=t2'; -dPjdTidVj=dPjdTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数 -dPdTidVj=dPidTidVj+dPjdTidVj;%最终非对角元素 @@ -%% deltaP/dVi_deltaThytaj 非对角元素 -t1=-ones(Busnum,1)*Volt.*real(GB).*(sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -dPidVidTj=t1-diag(diag(t1)); -t2=repmat(Init_Y',1,size(Init_Y,2)); -dPidVidTj=dPidVidTj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数 -t1=-ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat)); -dPjdVidTj=t1-diag(diag(t1)); -t2=t2'; -dPjdVidTj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数 -dPdVidTj=dPidVidTj+dPjdVidTj;%最终非对角元素 -%% deltaPi/dVi_deltaThyta 对角元素 -dPdVidTi=dPdTidVi;%最终对角元素 @@ -%% deltaP/dVi_dVj 非对角元素 -t1=-(real(GB).*(cos(AngleIJMat)+imag(GB).*sin(AngleIJMat))); -dPidVidVj=t1-diag(diag(t1));%去掉对角元素的 -t2=repmat(Init_Y',1,size(Init_Y,2)); -dPidVidVj=dPidVidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数 -t1=-(real(GB).*(cos(AngleIJMat)-imag(GB).*sin(AngleIJMat))); -dPjdVidVj=t1-diag(diag(t1)); -t2=t2'; -dPjdVidVj=dPjdVidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数 -dPdVidVj=dPidVidVj+dPjdVidVj;%最终非对角元素 @@ -%% deltaP/dVi_dVi 对角元素 -t1=diag(real(GB)); -t2=t1'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数 -dPidVidVi=-2*diag(t2); -dPidVjdVj=0; -dPdVidVi=dPidVidVi+dPidVjdVj;%最终对角元素 @@ -%% 生成APi -APi=zeros(2*Busnum,2*Busnum); -APi(1:2:2*Busnum,1:2:2*Busnum)=dPdTidTj;%%非对角 TT -APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVj;%%非对角 TV -APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTj;%%非对角 VT -APi(2:2:2*Busnum,2:2:2*Busnum)=dPdVidVj;%%非对角 VV -APi(1:2:2*Busnum,1:2:2*Busnum)=dPdTidTi;%%对角 -APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVi;%%对角 -APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTi;%%对角 -APi(2:2:2*Busnum,2:2:2*Busnum)=dPdVidVi;%%对角 -%% deltaQ/deltaThyta_deltaThyta 非对角元素 -t1=-Volt'*Volt; -t2=real(GB).*(sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -dQidTidTj=t1.*t2;%不去掉对角元素了，反正最后是要修正的 -t3=repmat(Init_Y',1,size(Init_Y,2)); -dQidTidTj=dQidTidTj.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数 -t2=real(GB).*(sin(AngleIJMat)+imag(GB).*cos(AngleIJMat)); -dQjdTidTj=t1.*t2; -t3=t3'; -dQjdTidTj=dQjdTidTj.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数 -dQdTidTj=dQidTidTj+dQjdTidTj;%最终非对角元素 -%% deltaQ/deltaThyta_deltaThyta 对角元素 -t1=dQidTidTj-diag(diag(dQidTidTj)); -t2=sum(t1,2); -t3=t2'.*Init_Y(2:2:size(Init_Y,2));%乘y的系数 -dQidTidTi=diag(t3); -t1=-Volt'*Volt; -t2=real(GB).*(sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -t3=t1.*t2; -t4=t3-diag(diag(t3)); -t5=sum(t4,1); -t6=t5.*Init_Y(2:2:size(Init_Y,2));%乘y的系数 -dQjdTidTi=diag(t6); -dQdTidTi=dQjdTidTi+dQidTidTi;%最终对角元素 -%% deltaQ/deltaThyta_deltaV 非对角元素 -t1=-Volt; -t2=real(GB).*(cos(AngleIJMat)+imag(GB).*sin(AngleIJMat)); -t3=ones(Busnum,1)*t1.*t2; -t4=repmat(Init_Y',1,size(Init_Y,2)); -dQidTidVj=t3.*t4(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数 -t2=real(GB).*(cos(AngleIJMat)-imag(GB).*sin(AngleIJMat)); -t3=Volt'*ones(1,Busnum).*t2; -t4=t4'; -dQjdTidVj=t2.*t3.*t4(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数 -dQdTidVj=dQidTidVj+dQjdTidVj;%最终非对角元素 -%% deltaQ/deltaThyta_deltaV 对角元素 -t1=sum(dQidTidVj,2)-diag(dQidTidVj); -t2=t1'.*Init_Y(2:2:size(Init_Y,2));%乘y的系数 -dQidTidVi=diag(t2); -t1=-Volt'*ones(1,Busnum).*(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat)); -t2=t1-diag(diag(t1)); -t3=sum(t2,1); -t4=t3.*Init_Y(2:2:size(Init_Y,2));%乘y的系数 -dQjdTidVi=diag(t4); -dQdTidVi=dQidTidVi+dQjdTidVi; -%% deltaQ/deltaV_deltaV 非对角元素 -t1=real(GB).*(sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -t2=repmat(Init_Y',1,size(Init_Y,2)); -t3=t1.*t2(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数 -dQidVidVj=t3; -t1=-real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -t2=t2'; -t3=t1.*t2(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数 -dQjdVidVj=t3; -dQdVidVj=dQidVidVj+dQjdVidVj; -%% deltaQ/deltaV_deltaV 对角元素 -t1=2*(real(GB).*(sin(AngleIJMat)-imag(GB).*cos(AngleIJMat))); -t2=diag(t1); -t3=t2'.*Init_Y(2:2:size(Init_Y,2));%乘y的系数 -dQidVidVi=diag(t3); -t1=-real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -t2=diag(t1); -t3=t1-diag(t2); -t4=sum(t3,1); -t5=t4.*Init_Y(2:2:size(Init_Y,2));%乘y的系数 -dQjdVidVi=diag(t5); -dQdVidVi=dQidVidVi+dQjdVidVi; -%% deltaQ/deltaV_deltaThyta 非对角元素 -t1=real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat); -t2=-ones(Busnum,1)*Volt.*(t1); -t3=repmat(Init_Y',1,size(Init_Y,2)); -t4=t2'.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数 -dQidVidTj=t4; -t1=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -t2=-ones(Busnum,1)*Volt.*(t1); -t3=t3'; -t4=t2.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数 -dQjdVidTj=t4; -dQdVidTj=dQidVidTj+dQjdVidTj; -%% deltaQ/deltaV_deltaThyta 对角元素 -dQdVidTi=dQdTidVi; -%% 生成AQi -AQi=zeros(2*Busnum,2*Busnum); -AQi(1:2:2*Busnum,1:2:2*Busnum)=dQdTidTj;%%非对角 TT -AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVj;%%非对角 TV -AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTj;%%非对角 VT -AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVj;%%非对角 VV -AQi(1:2:2*Busnum,1:2:2*Busnum)=dQdTidTi;%%对角 -AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVi;%%对角 -AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTi;%%对角 -AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVi;%%对角 -%% 生成ddf -t=[zeros(2*14-2*5,2*14); - zeros(2*5,2*14-2*5),AQi+APi; -]; -ddf=t; -end \ No newline at end of file diff --git a/func_ddg.asv b/func_ddg.asv deleted file mode 100644 index a2e74a5..0000000 --- a/func_ddg.asv +++ /dev/null @@ -1,14 +0,0 @@ -function ddg=func_ddg(PGi,PVi,Busnum,ContrlCount,RestraintCount,Loadi,PD,QD) -% 计算系数 -InitZW=Init_Z+InitW; -loc=InitZW(length(PGi)+length(PVi)+length(Loadi)*2:RestraintCount); -PD_=PD(Loadi); -QD_=QD(Loadi) -t1=sparse(length(PGi)+length(PVi),ContrlCount); -t21=sparse(length(PD_),length(PGi)+length(PVi)); -t22=sparse(-2./PD_.^2-2*(QD_.^2-PD_.^2)./( PD_.^2+QD_.^2).^2)*diag(loc);%dF/dPdP -t23=sparse(4*PD_.*PQ_./(PD_.^2+PQ_.^2).^2); %dF/dPdQ -t24=sparse(length(PD_),2*Busnum); -t31=t21; -t32=sparse(-2.) -end \ No newline at end of file diff --git a/func_ddg.m b/func_ddg.m deleted file mode 100644 index ead0e3c..0000000 --- a/func_ddg.m +++ /dev/null @@ -1,23 +0,0 @@ -function ddg=func_ddg(PGi,PVi,Busnum,ContrlCount,RestraintCount,Loadi,PD,QD,Init_Z,Init_W) -% 计算系数 -InitZW=Init_Z+Init_W; -loc=InitZW(length(PGi)+length(PVi)+length(Loadi)*1+Busnum+1:RestraintCount); -PD_=PD(Loadi); -QD_=QD(Loadi); -t1=sparse(length(PGi)+length(PVi),ContrlCount); -t21=sparse(length(PD_),length(PGi)+length(PVi)); -t22=diag(sparse(-2./PD_.^2-2*(QD_.^2-PD_.^2)./( PD_.^2+QD_.^2).^2))*diag(loc);%dF/dPdP -t23=diag(sparse(4*PD_.*QD_./(PD_.^2+QD_.^2).^2))*diag(loc); %dF/dPdQ -t24=sparse(length(PD_),2*Busnum); -t31=t21; -t32=t23; -t33=diag(sparse(-2*(PD_.^2-QD_.^2)./(PD_.^2+QD_.^2).^2 ))*diag(loc); -t34=t24; -t4=sparse(Busnum*2,ContrlCount); -ddg=[ - t1; - t21,t22,t23,t24; - t31,t32,t33,t34; - t4; -]; -end \ No newline at end of file diff --git a/func_ddh.m b/func_ddh.m deleted file mode 100644 index 8403554..0000000 --- a/func_ddh.m +++ /dev/null @@ -1,58 +0,0 @@ -function ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount) -%决定用循环重写 -%ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2; - -%AngleIJ=AngleIJMat-angle(GB); -mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); -mat_INV_AngleIJ=mat_AngleIJ'; -yP=Init_Y(1:size(Init_Y,2)/2);%暂时改这里 20111227 -yQ=Init_Y(size(Init_Y,2)/2+1:size(Init_Y,2));%暂时改这里 20111227 -t1=-sparse(1:Busnum,1:Busnum,Y.*cos(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP',Busnum,Busnum); -t2=sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP',Busnum,Busnum)*Y.*cos(mat_AngleIJ); -t3=(t1+t2)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); -t4=-(sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt',Busnum,Busnum) -sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP',Busnum,Busnum); -ddPdTdT=t3+t4;%ok1 -t1=(-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt',Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,yP,Busnum,Busnum); -t2= -sparse(1:Busnum,1:Busnum, sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP' ,Busnum,Busnum)*Y.*sin(mat_AngleIJ)+sparse(1:Busnum,1:Busnum,Y.*sin(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP',Busnum,Busnum); -ddPdVdT=t1+t2;%ok1 -t1=diag( Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yP'); -t2=diag(yP)*Y.*sin(mat_AngleIJ)*diag(Volt); -t3=-diag(yP)*diag(Y.*sin(mat_AngleIJ)*Volt'); -t4=-Y.*sin(mat_INV_AngleIJ)*diag( diag(Volt)*yP' ); -ddPdTdV=t1+t2+t3+t4;%存疑与我的不一样 -t1=Y.*cos(mat_INV_AngleIJ)*diag(yP); -t2=diag(yP)*Y.*cos(mat_AngleIJ); -ddPdVdV=t1+t2; -t1=-diag(Y.*sin(mat_AngleIJ)*Volt'); -t2=diag(Volt)*Y.*sin(mat_INV_AngleIJ); -t3=(t1+t2)*diag( diag(Volt)*yQ' ); -t4=-diag( diag(Volt)*yQ' )*Y.*sin(mat_AngleIJ); - -t5=diag(Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yQ'); -t6=-(t4+t5)*diag(Volt); -ddQdTdT=t3+t6;%ok1 -t1=(diag(Y.*cos(mat_AngleIJ)*Volt')-diag(Volt)*Y.*cos(mat_INV_AngleIJ) )*diag(yQ); -t2=+diag( diag(Volt)*yQ' )*Y.*cos(mat_AngleIJ)-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ'); -ddQdVdT=t1+t2; -t1=Y.*cos(mat_INV_AngleIJ)*diag(diag(Volt)*yQ'); -t2=diag(yQ)*diag(Y.*cos(mat_AngleIJ)*Volt'); -t3=-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ'); -t4=-diag(yQ)*Y.*cos(mat_AngleIJ)*diag(Volt); -ddQdTdV=t1+t2+t3+t4; -t1=Y.*sin(mat_INV_AngleIJ)*diag(yQ); -t2=diag(yQ)*Y.*sin(mat_AngleIJ); -ddQdVdV=t1+t2; -% t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV; -% ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT; -% ]; -t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT; - ddPdTdV+ddQdTdV,ddPdTdT+ddQdTdT; -]; -sizePGi=size(PGi,1); -sizePVi=size(PVi,1); -sizeLoadi=size(Loadi,1)*2; -ddh=[ - sparse(sizePGi+sizePVi+sizeLoadi,ContrlCount); - sparse(2*Busnum,sizePVi+sizePGi+sizeLoadi),-t; - ]; -end \ No newline at end of file diff --git a/func_deltF.asv b/func_deltF.asv deleted file mode 100644 index 8fb2e42..0000000 --- a/func_deltF.asv +++ /dev/null @@ -1,19 +0,0 @@ -function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum) -%t1=PG(setdiff(PVi,Balance)); -% t2=Volt'*Volt; -% t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -% t4=t2.*t3; -% t5=sum(t4,2); -% PBal=t5(Balance); -% PPG=([PQ(1),PBal])';%暂时用土办法处理一下 -%% -t1=2 -t2=2*wD.*(PD-PD0); -deltF=[ - zeros(size(PGi)); - zeros(size(PVi)); - t2; - zeros(2*Busnum,1); -]; - -end \ No newline at end of file diff --git a/func_deltF.m b/func_deltF.m deleted file mode 100644 index 2212075..0000000 --- a/func_deltF.m +++ /dev/null @@ -1,15 +0,0 @@ -function deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD0,QD,Busnum,Loadi) - -t1=2*wPG.*(PG(PGi) - PG0(PGi) ); -t2=2*wQG.*(QG(PVi) - QG0(PVi) ); -t3=2*wPD.*(PD(Loadi)-PD0(Loadi)); -t4=2*wQD.*(QD(Loadi)-QD0(Loadi)); -deltF=[ - sparse(t1); - sparse(t2); - sparse(t3); - sparse(length(Loadi),1); - sparse(2*Busnum,1); -]; - -end \ No newline at end of file diff --git a/func_deltG.asv b/func_deltG.asv deleted file mode 100644 index 7048616..0000000 --- a/func_deltG.asv +++ /dev/null @@ -1,52 +0,0 @@ -function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD) -%% -PD_=PD(Loadi); -QD_=QD(Loadi); -sizePGi=size(PGi,1); -sizePVi=size(PVi,1); -sizeLoadi=size(Loadi,1); -%% -dg1_dPg=sparse(1:sizePGi,1:sizePGi,ones(sizePGi,1),sizePGi,sizePGi); -dg2_dPg=sparse(sizePGi,sizePVi); -dg3_dPg=sparse(sizePGi,sizeLoadi); -dg4_dPg=sparse(sizePGi,length(Loadi)); -dg5_dPg=sparse(sizePGi,Busnum); -dg6_dPg=sparse(sizePGi,length(Loadi)); -%% -dg1_dQr=sparse(sizePVi,sizePGi); -dg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi); -dg3_dQr=sparse(sizePVi,sizeLoadi); -dg4_dQr=sparse(sizePVi,length(Loadi)); -dg5_dQr=sparse(sizePVi,Busnum); -dg6_dQr=sparse(sizePVi,length(Loadi)); -%% -dg1_dPD=sparse(size(Loadi,1),size(PGi,1)); -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),length(Loadi)); -dg5_dPD=sparse(size(Loadi,1),Busnum); -dg6_dPD=diag(2./PD_-2*PD_./(PD_.^2+QD_.^2)); -%% -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(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1)); -dg5_dQD=sparse(size(Loadi,1),Busnum); -dg6_dQD=diag(-2*QD_./(PD_.^2+QD_.^2)); -%% -dg1_dx=sparse(2*Busnum,sizePGi); -dg2_dx=sparse(2*Busnum,sizePVi); -dg3_dx=sparse(2*Busnum,sizeLoadi); -dg4_dx=sparse(2*Busnum,length(Loadi)); -dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum); - sparse(Busnum,Busnum); - ]; -dg6_dx=sparse(2*Busnum,Busnum); -%% -deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg4_dPg,dg5_dPg,dg6_dPg; - dg1_dQr,dg2_dQr,dg3_dQr,dg4_dQr,dg5_dQr,dg6_dQr; - dg1_dPD,dg2_dPD,dg3_dPD,dg4_dPD,dg5_dPD,dg6_dPD; - dg1_dQD,dg2_dQD,dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD; - dg1_dx,dg2_dx,dg3_dx,dg4_dx,dg5_dx,dg6_dx; -]; -end \ No newline at end of file diff --git a/func_deltG.m b/func_deltG.m deleted file mode 100644 index 26fc17a..0000000 --- a/func_deltG.m +++ /dev/null @@ -1,58 +0,0 @@ -function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD) -%% -PD_=PD(Loadi); -QD_=QD(Loadi); -sizePGi=size(PGi,1); -sizePVi=size(PVi,1); -sizeLoadi=size(Loadi,1); -%% -dg1_dPg=sparse(1:sizePGi,1:sizePGi,ones(sizePGi,1),sizePGi,sizePGi); -dg2_dPg=sparse(sizePGi,sizePVi); -dg3_dPg=sparse(sizePGi,sizeLoadi); -%dg4_dPg=sparse(sizePGi,length(Loadi)); -dg5_dPg=sparse(sizePGi,Busnum); -dg6_dPg=sparse(sizePGi,length(Loadi)); -%% -dg1_dQr=sparse(sizePVi,sizePGi); -dg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi); -dg3_dQr=sparse(sizePVi,sizeLoadi); -%dg4_dQr=sparse(sizePVi,length(Loadi)); -dg5_dQr=sparse(sizePVi,Busnum); -dg6_dQr=sparse(sizePVi,length(Loadi)); -%% -dg1_dPD=sparse(size(Loadi,1),size(PGi,1)); -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),length(Loadi)); -dg5_dPD=sparse(size(Loadi,1),Busnum); -dg6_dPD=diag(2./PD_-2*PD_./(PD_.^2+QD_.^2)); -%% -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(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1)); -dg5_dQD=sparse(size(Loadi,1),Busnum); -dg6_dQD=diag(-2*QD_./(PD_.^2+QD_.^2)); -%% -dg1_dx=sparse(2*Busnum,sizePGi); -dg2_dx=sparse(2*Busnum,sizePVi); -dg3_dx=sparse(2*Busnum,sizeLoadi); -%dg4_dx=sparse(2*Busnum,length(Loadi)); -dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum); - sparse(Busnum,Busnum); - ]; -dg6_dx=sparse(2*Busnum,length(Loadi)); -%% -% deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg4_dPg,dg5_dPg,dg6_dPg; -% dg1_dQr,dg2_dQr,dg3_dQr,dg4_dQr,dg5_dQr,dg6_dQr; -% dg1_dPD,dg2_dPD,dg3_dPD,dg4_dPD,dg5_dPD,dg6_dPD; -% dg1_dQD,dg2_dQD,dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD; -% dg1_dx,dg2_dx,dg3_dx,dg4_dx,dg5_dx,dg6_dx; -% ]; -deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg5_dPg,dg6_dPg; - dg1_dQr,dg2_dQr,dg3_dQr,dg5_dQr,dg6_dQr; - dg1_dPD,dg2_dPD,dg3_dPD,dg5_dPD,dg6_dPD; - dg1_dQD,dg2_dQD,dg3_dQD,dg5_dQD,dg6_dQD; - dg1_dx,dg2_dx,dg3_dx,dg5_dx,dg6_dx; -]; -end \ No newline at end of file diff --git a/func_deltH.asv b/func_deltH.asv deleted file mode 100644 index 317b2f3..0000000 --- a/func_deltH.asv +++ /dev/null @@ -1,8 +0,0 @@ -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_dx = jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c); %形成雅克比矩阵 -deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx']; -end \ No newline at end of file diff --git a/func_deltH.m b/func_deltH.m deleted file mode 100644 index 3c0fb7e..0000000 --- a/func_deltH.m +++ /dev/null @@ -1,10 +0,0 @@ -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_dQD;dH_dx']; - -end \ No newline at end of file diff --git a/func_deltdeltF.asv b/func_deltdeltF.asv deleted file mode 100644 index 87f3a1a..0000000 --- a/func_deltdeltF.asv +++ /dev/null @@ -1,11 +0,0 @@ -function deltdeltF=func_deltdeltF(Busnum,PVi,PGi,wG,wD,ContrlCount) - -ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数 -C=[wG' zeros(size(PVi))' wD']; -sizeC=size(C,2); -diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC); -deltdeltF=[ - diagC*2,sparse(sizeC,ContrlCount-sizeC); - sparse(ContrlCount-sizeC,ContrlCount); - ]; -end \ No newline at end of file diff --git a/func_deltdeltF.m b/func_deltdeltF.m deleted file mode 100644 index 820fc12..0000000 --- a/func_deltdeltF.m +++ /dev/null @@ -1,11 +0,0 @@ -function deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount) - -%ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数 -C=[wPG' wQG' wPD' wQD']; -sizeC=size(C,2); -diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC); -deltdeltF=[ - diagC*2,sparse(sizeC,ContrlCount-sizeC); - sparse(ContrlCount-sizeC,ContrlCount); - ]; -end \ No newline at end of file diff --git a/jacobian_M.m b/jacobian_M.m deleted file mode 100644 index c224cf2..0000000 --- a/jacobian_M.m +++ /dev/null @@ -1,20 +0,0 @@ -function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c) -%************************************************************************** -% 程序功能 : 子函数——形成雅可比矩阵Jacobian -% 编 者： -% 编制时间：2010.12 -%************************************************************************** -%%参照图书馆6楼的书编写 -%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q -AngleIJ=UAngel(r)-UAngel(c)-Angle'; -mat_AngleIJ=sparse(r,c,AngleIJ,Busnum,Busnum); -mat_IvAngleIJ=mat_AngleIJ'; -H=sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt',Busnum,Busnum)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); -N=-sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt',Busnum,Busnum)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); -J=sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt',Busnum,Busnum)+Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); -L=sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt',Busnum,Busnum)+Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); -t1=[J,L; - H,N; -]'; -Jacob=-t1; -end \ No newline at end of file diff --git a/modifyadmmatrix.asv b/modifyadmmatrix.asv deleted file mode 100644 index a917e8d..0000000 --- a/modifyadmmatrix.asv +++ /dev/null @@ -1,18 +0,0 @@ -function [new_G,new_B,GB,Y,r,c,Angle] = modifyadmmatrix(ii,jj,G,GG,B,BB) -%************************************************************************** -% 程序功能 : 子函数——形成节点导纳矩阵Y -% 编 者： -% 编制时间：2010.12 -%************************************************************************** -%% 支路导纳计算 -new_G=G; -new_G(ii,jj)=new_G(ii,jj)-G(ii,jj); -new_G(jj,ii)=new_G(jj,ii)-G(jj,ii); -new_G(ii,ii)=new_G(ii,ii)+G(ii,jj); -new_G(jj,jj)=new_G(jj,jj)+G(ii,jj); - -%% 化作极坐标形式 -GB = G+B.*1i; %将电导，电纳合并，写成复数形式 -Y = abs(GB); %求节点导纳幅值 -[r,c] = find(Y); -Angle = angle(GB(GB~=0)); %求节点导纳角度 \ No newline at end of file diff --git a/modifyadmmatrix.m b/modifyadmmatrix.m deleted file mode 100644 index de48607..0000000 --- a/modifyadmmatrix.m +++ /dev/null @@ -1,23 +0,0 @@ -function [new_G,new_B,GB,Y,r,c,Angle] = modifyadmmatrix(ii,jj,G,B) -%************************************************************************** -% 程序功能 : 子函数——形成节点导纳矩阵Y -% 编 者： -% 编制时间：2010.12 -%************************************************************************** -%% 支路导纳计算 -new_G=G; -new_G(ii,jj)=new_G(ii,jj)-G(ii,jj); -new_G(jj,ii)=new_G(jj,ii)-G(jj,ii); -new_G(ii,ii)=new_G(ii,ii)+G(ii,jj); -new_G(jj,jj)=new_G(jj,jj)+G(ii,jj); -new_B=B; -new_B(ii,jj)=new_B(ii,jj)-B(ii,jj); -new_B(jj,ii)=new_B(jj,ii)-B(jj,ii); -new_B(ii,ii)=new_B(ii,ii)+B(ii,jj); -new_B(jj,jj)=new_B(jj,jj)+B(ii,jj); - -%% 化作极坐标形式 -GB = new_G+new_B.*1i; %将电导，电纳合并，写成复数形式 -Y = abs(GB); %求节点导纳幅值 -[r,c] = find(Y); -Angle = angle(GB(GB~=0)); %求节点导纳角度 \ No newline at end of file diff --git a/pf.asv b/pf.asv deleted file mode 100644 index 575e701..0000000 --- a/pf.asv +++ /dev/null @@ -1,34 +0,0 @@ -function [kmax,Precision,Uangle,U,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(FileName) -%************************************************************************** -% 程序名称：电力系统潮流计算程序 -% 程序算法：极坐标下的牛顿-拉夫逊法 -% 程序功能：主函数 -% 程序编者： -% 编制时间：2010.12 -%************************************************************************** -clc; -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]= openfile(FileName); -%% 形成节点导纳矩阵 -[GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,... - Transforx,Transfork0,Branchi,Branchb); -[P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum); %求功率不平衡量 -disp('迭代次数i 最大不平衡量'); -%% 循环体计算 -for i = 0:kmax - [Jacob,PQ,U,Uangle] = jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c); %形成雅克比矩阵 - % disp('第一次雅克比'); - m = max(abs(PQ)); - m=full(m); - fprintf(' %u %.8f \n',i,m); - if m > Precision %判断不平衡量是否满足精度要求 - [Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U); %求解修正方程，更新电压变量 - else - disp(['收敛，迭代次数为',num2str(i),'次']); - break %若满足精度要求，则计算收敛 - end -end -toc; -end \ No newline at end of file diff --git a/test.m b/test.m deleted file mode 100644 index e7b3690..0000000 --- a/test.m +++ /dev/null @@ -1,97 +0,0 @@ -%% Test -clear -arraymaxPDError=zeros(54,1); -for HH=0:53 - - arraymaxPDError(HH+1)=-100; - for CC=1:250 - tic - [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('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); - %pf('c:/file31.txt'); - %pf('ieee10471PG.dat'); - - %PVi电压节点序号 - %PVu电压节点电压标幺值 - Volt; - UAngel*180/3.1415926; - %Precision=Precision/1000000; - %Precision=Precision*1000; - %% 通过潮流计算PG - AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); - PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt'; - - %% 初值 - PG0=PG; - PD0=PD; - %% - PG0(Balance)=PGBal(Balance); - %% - [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,HH); - 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; - 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; - end - plotGap(KK+1)=Gap; - Init_u=Gap/2/RestraintCount*CenterA; - %AngleIJMat=0; - %% 开始计算OPF - %% 形成等式约束的雅克比 - deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi); - %% 形成不等式约束的雅克比 - deltG=func_deltG(Busnum,PVi,PGi,Loadi); - %% - L_1Z=diag(Init_Z./Init_L); - U_1W=diag(Init_W./Init_U); - %% 形成海森阵 - deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount); - %% 形成ddHy - ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount); - %% 开始构建ddg - ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi); - %% 开始构建deltF - deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi); - - %% 形成方程矩阵 - 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,Loadi); - Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND); - Ly=Mat_H; - Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK); - Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); - %%YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); - %% 开始解方程 - 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); - Gap=(Init_L*Init_Z'-Init_U*Init_W'); - KK=KK+1; - end - fprintf('迭代次数%d\n',KK); - ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi); - %DrawGap(plotGap); - %% - %Volt=full(Volt'); - %PD=full(PD); - %% 统计PD误差 - maxPDError=max(abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) )); - if maxPDError>arraymaxPDError(HH+1) - arraymaxPDError(HH+1)=maxPDError; - end - toc - end -end - -