diff --git a/0PD0.mat b/0PD0.mat deleted file mode 100644 index 00bfc6c..0000000 Binary files a/0PD0.mat and /dev/null differ diff --git a/0QD0.mat b/0QD0.mat deleted file mode 100644 index 733cc73..0000000 Binary files a/0QD0.mat and /dev/null differ diff --git a/12PD0.mat b/12PD0.mat deleted file mode 100644 index 65f1f29..0000000 Binary files a/12PD0.mat and /dev/null differ diff --git a/12QD0.mat b/12QD0.mat deleted file mode 100644 index 371b4f2..0000000 Binary files a/12QD0.mat and /dev/null differ diff --git a/16PD0.mat b/16PD0.mat deleted file mode 100644 index e4e2f65..0000000 Binary files a/16PD0.mat and /dev/null differ diff --git a/16QD0.mat b/16QD0.mat deleted file mode 100644 index fe41c6b..0000000 Binary files a/16QD0.mat and /dev/null differ diff --git a/20PD0.mat b/20PD0.mat deleted file mode 100644 index f8e0b19..0000000 Binary files a/20PD0.mat and /dev/null differ diff --git a/20QD0.mat b/20QD0.mat deleted file mode 100644 index 0296b5f..0000000 Binary files a/20QD0.mat and /dev/null differ diff --git a/4PD0.mat b/4PD0.mat deleted file mode 100644 index 3c6017a..0000000 Binary files a/4PD0.mat and /dev/null differ diff --git a/4QD0.mat b/4QD0.mat deleted file mode 100644 index 5d5ef3a..0000000 Binary files a/4QD0.mat and /dev/null differ diff --git a/8PD0.mat b/8PD0.mat deleted file mode 100644 index 5d83ebe..0000000 Binary files a/8PD0.mat and /dev/null differ diff --git a/8QD0.mat b/8QD0.mat deleted file mode 100644 index c73bbab..0000000 Binary files a/8QD0.mat and /dev/null differ diff --git a/@ForThesis/MaxBranchDeviation.m b/@ForThesis/MaxBranchDeviation.m index ac15e4c..b715597 100644 --- a/@ForThesis/MaxBranchDeviation.m +++ b/@ForThesis/MaxBranchDeviation.m @@ -29,8 +29,8 @@ end Volt0=U; UAngel0=Uangle; -[dispLineloss0 dispTransloss0]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0); -[dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt,UAngel); +[dispLineloss0 dispTransloss0]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0);%测量值 +[dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt,UAngel);%估计值 t1=(dispLineloss0(:,3) - dispLineloss(:,3))./dispLineloss0(:,3); t2=(dispTransloss0(:,3) - dispTransloss(:,3))./dispTransloss0(:,3); t11=dispLineloss0(:,3)>1e-5;% 太小的值不计算 diff --git a/FormG.m b/FormG.m index f137773..7cf7e72 100644 --- a/FormG.m +++ b/FormG.m @@ -1,8 +1,10 @@ -function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi) +function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi,Vbi) Mat_G=[ sparse(PD(Loadi)); sparse(QD(Loadi)); Volt'; + Volt'; + Vbi; ]; end \ No newline at end of file diff --git a/FormH.m b/FormH.m index 0779496..e655496 100644 --- a/FormH.m +++ b/FormH.m @@ -1,10 +1,4 @@ -function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND,Loadi) -%% -%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)); +function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi) AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt'; dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt'; diff --git a/FormLw.asv b/FormLw.asv deleted file mode 100644 index 987b6ee..0000000 --- a/FormLw.asv +++ /dev/null @@ -1,25 +0,0 @@ -function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity) -KK=999; -%PU=GenU(:,2);%发电机有功上界 -PU=5*ones(length(GenU(:,2)),1); -%QU=PVQU(:,1);%发电机无功上界 -QU=5*ones(length(PVQU(:,1)),1); -VoltU=(1.1)*ones(1,Busnum); -%VoltU=10*ones(1,Busnum); -% PDU=PD0(Loadi); -PDU=noDataTransCapacity; -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=0.85; -% QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF; -t1=([PU',QU',PDU',QDU',VoltU])'; -t2=Mat_G+Init_U'-t1; -Lw=t2; - -end \ No newline at end of file diff --git a/FormLw.m b/FormLw.m index 8c02646..1a7b62e 100644 --- a/FormLw.m +++ b/FormLw.m @@ -9,14 +9,22 @@ PDU=PD0(Loadi); PDU(PDU>0)=1.200*PDU(PDU>0); PDU(PDU<0)=0.800*PDU(PDU<0); PDU(PDU==0)=0.400; +realPD=PD0(Loadi); +indPD=find(realPD>0); +PDU(indPD(3:12:end))=1.55*realPD(indPD(3:12:end)); +PDU(indPD(9:12:end))=1.05*realPD(indPD(9:12:end)); %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.400; +realQD=QD0(Loadi); +indQD=find(realQD>0); +QDU(indQD(3:12:end))=1.55*realQD(indQD(3:12:end)); +QDU(indQD(9:12:end))=1.05*realQD(indQD(9:12:end)); % PF=0.85; % QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF; -t1=([PDU',QDU',VoltU])'; +t1=([PDU',QDU',VoltU,VoltU,0*ones(1,Busnum)])'; t2=Mat_G+Init_U'-t1; Lw=t2; diff --git a/FormLz.asv b/FormLz.asv deleted file mode 100644 index 48a736e..0000000 --- a/FormLz.asv +++ /dev/null @@ -1,20 +0,0 @@ -function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF) -KK=999; - -VoltL=(0.9)*ones(1,Busnum); -%VoltL=-10*ones(1,Busnum); -PDL=PD0(Loadi); -PDL(PDL>0)=0.700*PDL(PDL>0); -PDL(PDL<0)=1.300*PDL(PDL<0); -PDL(PDL==0)=-0.400; -%PDL=-10*ones(length(Loadi),1); -QDL=QD0(Loadi); -QDL(QDL>0)=0.700*QDL(QDL>0); -QDL(QDL<0)=1.300*QDL(QDL<0); -QDL(QDL==0)=-0.400; -% QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF; -t1=([PDL',QDL',VoltL])'; -t2=Mat_G-Init_L'-t1; -Lz=t2; - -end \ No newline at end of file diff --git a/FormLz.m b/FormLz.m index 1f58ebc..b2b8984 100644 --- a/FormLz.m +++ b/FormLz.m @@ -7,13 +7,21 @@ PDL=PD0(Loadi); PDL(PDL>0)=0.800*PDL(PDL>0); PDL(PDL<0)=1.200*PDL(PDL<0); PDL(PDL==0)=-0.400; +realPD=PD0(Loadi); +indPD=find(realPD>0); +PDL(indPD(3:12:end))=0.95*realPD(indPD(3:12:end)); +PDL(indPD(9:12:end))=0.45*realPD(indPD(9:12:end)); %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.400; +realQD=QD0(Loadi); +indQD=find(realQD>0); +QDL(indQD(3:12:end))=0.95*realQD(indQD(3:12:end)); +QDL(indQD(9:12:end))=0.95*realQD(indQD(9:12:end)); % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF; -t1=([PDL',QDL',VoltL])'; +t1=([PDL',QDL',VoltL,VoltL,0*ones(1,Busnum)])'; t2=Mat_G-Init_L'-t1; Lz=t2; diff --git a/JSMJZM.asv b/JSMJZM.asv deleted file mode 100644 index 86be91d..0000000 --- a/JSMJZM.asv +++ /dev/null @@ -1,12 +0,0 @@ -JSM=[0.0575 0.051 0.0029]; -JZM=[0.003 0.02629 0.0582]; -plot(1:3,JSM,'k'); -axis([0 4 0 0.07]) -text(1,JSM(1)+0.001,'0.0575') -text(2,JSM(2)+0.001,'0.051') -text(3,JSM(3)+0.001,'0.0029') -hold on -plot(1:3,JZM,'k--'); -text(1,JZM(1)+0.001,'0.003') -text(2,JZM(2)+0.001,'0.02629') -text(3,JZM(3)+0.001,'0.0029') \ No newline at end of file diff --git a/JSMJZM.m b/JSMJZM.m index 936502b..ce51633 100644 --- a/JSMJZM.m +++ b/JSMJZM.m @@ -1,18 +1,18 @@ -JSM=[0.0575 0.051 0.0029]; -JZM=[0.0003 0.02629 0.0582]; -for I=0.01:0.01:0.06 +JSM=[0.070244 0.0533 0.004379]; +JZM=[0.002275 0.047311 0.070318]; +for I=0.01:0.01:0.08 line([0.5 3.5],[I I],'Color',[220 220 220]/255) end hold on plot(1:3,JSM,'k'); ylabel('统计误差') -xlabel('情况') -axis([0.5 3.5 0 0.065]) -text(1,JSM(1)+0.001,'0.0575') -text(2,JSM(2)+0.001,'0.0510') -text(3,JSM(3)+0.001,'0.0029') +xlabel('情形') +axis([0.5 3.5 0 0.075]) +text(1,JSM(1)+0.001,'0.0702') +text(2,JSM(2)+0.001,'0.0533') +text(3,JSM(3)+0.001,'0.0043') hold on plot(1:3,JZM,'k--'); -text(1-.3,JZM(1)+0.001,'0.0003') -text(2,JZM(2)-0.001,'0.0263') -text(3,JZM(3)-0.001,'0.0029') +text(1-.3,JZM(1)+0.001,'0.0023') +text(2-0.3,JZM(2)-0.001,'0.0473') +text(3,JZM(3)-0.001,'0.0703') diff --git a/JSNZN.asv b/JSNZN.asv deleted file mode 100644 index ea4d42c..0000000 --- a/JSNZN.asv +++ /dev/null @@ -1,14 +0,0 @@ -JSN=[0.1213 0.1092 0.0587]; -JZN=[0.0971 0.1314 0.1408]; -plot(1:3,JSM,'k'); -ylabel('值') -xlabel('情况') -axis([0.5 3.5 0 0.065]) -text(1,JSM(1)+0.001,'0.1213') -text(2,JSM(2)+0.001,'0.0510') -text(3,JSM(3)+0.001,'0.0029') -hold on -plot(1:3,JZM,'k--'); -text(1-.3,JZM(1)+0.001,'0.0003') -text(2,JZM(2)-0.001,'0.0263') -text(3,JZM(3)-0.001,'0.0029') \ No newline at end of file diff --git a/Lineloss.asv b/Lineloss.asv deleted file mode 100644 index 16e5eaa..0000000 --- a/Lineloss.asv +++ /dev/null @@ -1,38 +0,0 @@ -%% 计算线损 -function [dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,Branchi,Branchg,Branchb,k0,Volt,Angle) -%format long -% fprintf('功率为有名值\n'); -% fprintf('节点号\t节点号\t有功损耗 MW\t无功损耗 MVar') -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; - -%% 另一种计算方式begin -ss=1*(Volt(Linei)'.^2.*abs(yij).*cos( angle(yij) ) -Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)' - angle(yij)).*abs(yij)); -ss=(Volt(Linei)'.^2+Volt(Linej)'.^2).*abs(yij).*cos(angle(yij))-2*Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)').*cos( - angle(yij)).*abs(yij); -ss=Volt(Linei)'.*Volt(Linej)'.*abs(yij).*cos() -%% 另一种计算方式end -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; -%% 接地支路损耗 -% 没有考虑变压器变比 -deltTransS =deltTransS+sum(cmpVolt(Branchi).*conj((cmpVolt(Branchi).*(Branchg+1j*Branchb)))); -%% -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 index 615244e..b1d91e5 100644 --- a/Modification.m +++ b/Modification.m @@ -1,4 +1,4 @@ -function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,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) +function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi]=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,Vbi) AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU); %fprintf('AlphaP %f\n',full(AlphaP)); AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW); @@ -24,4 +24,5 @@ balVolt=Volt(Balance); Volt=Volt+AlphaP*t(1:Busnum); Volt(Balance)=balVolt; UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum); +Vbi=Vbi+AlphaP*t(2*Busnum+1:end)'; end \ No newline at end of file diff --git a/OPF.asv b/OPF.asv deleted file mode 100644 index fe989df..0000000 --- a/OPF.asv +++ /dev/null @@ -1,137 +0,0 @@ -tic -clc -clear -thesis=ForThesis(4,8); -for II=1:4 - [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('c:/newFIle.txt'); - %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; - %% - PG0(Balance)=PGBal(Balance); - PG(Balance)=PGBal(Balance); - QG0(Balance)=QGBal(Balance); - QG0(PVi)=QGBal(PVi); - QG(PVi)=QGBal(PVi); - %% 真实值 - RealPG=PG0; - RealQG=QG0; - RealPD=PD0; - RealQD=QD0; - %% - %PGi=zeros(1,1); - [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,RealPD,RealQD,QD,PD); - Gap=(Init_L*Init_Z'-Init_U*Init_W'); - KK=0; - plotGap=zeros(1,60); - ContrlCount=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); - %% - Precision=Precision/10; - - %% 加误差 - PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1)); - QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1)); -% PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1)); -% QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1)); - %% 读变压器容量 - %[noDataTransNum noDataTransCapacity noDataTransPowerFactor]=ReadNoDataTrans('C:/b/东际911_2751267_2012-09-05/iPso_东际911_2751267_2012-09-05_变压器无负载.txt'); - noDataTransCapacity=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,RestraintCount,Loadi,PD,QD); - %% 开始构建deltF - deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,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,PD,RealPD,RealQD,Loadi,KK,PF); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF,noDataTransCapacity); - 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)-RealPD(Loadi))./RealPD(Loadi) ); - maxPDError=max(absPDLoad(absPDLoad<10)) - Loadi(maxPDError==absPDLoad) - absQDLoad=abs( (QD(Loadi)-RealQD(Loadi))./RealQD(Loadi) ); - maxQDError=max(absQDLoad(absQDLoad<10)) - Loadi(maxQDError==absQDLoad) - disp('index'); - %Loadi(absPDLoad==maxPDError); - %% 计算总线损 - totalLoss=(sum(PG)-sum(PD(Loadi)))*100; - fprintf('总的损耗为%f(MW 有名值)\n',full(totalLoss)); - fprintf('线损率为 %f\n',full(totalLoss/sum(PG))); - %% 计算各线损 - %Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel); - thesis=thesis.AddPDQDPGQG(PD(Loadi),QD(Loadi),PG(Balance),QG(PVi)); - -end -PD(Loadi)=thesis.MeanPD(); -QD(Loadi)=thesis.MeanQD(); -PG(Balance)=thesis.MeanPG(); -QG(PVi)=thesis.MeanQG(); -thesis.MaxDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi)) -thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi)); -toc - diff --git a/OPF.m b/OPF.m index 5a67bfc..162384d 100644 --- a/OPF.m +++ b/OPF.m @@ -3,146 +3,95 @@ clc clear %% 存在问题 % 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123 -%% -thesis=ForThesis(1,45); -for II=1:1 - [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,Branchi,Branchg,Branchb,Transfork0]= ... - pf('E:/算例/柳金Ⅰ926_21671693_2012-09-06/newFIle20 - 计算用2 - (最终用这个做20点的计算).txt'); +%% +thesis=ForThesis(1,62); +[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,Branchi,Branchg,Branchb,Transfork0]= ... + pf('C:\lz\结题后重新计算\重构结果\78\木材厂堤Ⅰ线_2229880_2013-10-29\木材厂堤Ⅰ线_2229880_2013-10-29_iPso_newFile.txt'); % pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.txt'); - %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; - Volt0=Volt; - UAngel0=UAngel; - %% - PG0(Balance)=PGBal(Balance); - PG(Balance)=PGBal(Balance); - QG0(Balance)=QGBal(Balance); - QG0(PVi)=QGBal(PVi); - QG(PVi)=QGBal(PVi); - %% 真实值 - RealPG=PG0; - RealQG=QG0; - RealPD=PD0; - RealQD=QD0; - %% - %PGi=zeros(1,1); - [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,RealPD,RealQD,QD,PD); - Gap=(Init_L*Init_Z'-Init_U*Init_W'); - KK=0; - plotGap=zeros(1,60); - ContrlCount=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); - %% - Precision=Precision/1; - %% 加误差 - PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1)); - QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1)); -% save('20PD0.mat','PD0'); -% save('20QD0.mat','QD0'); - load('20PD0.mat'); - load('20QD0.mat'); -% PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1)); -% QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1)); - %% 读变压器容量 - noDataTransCapacity=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,RestraintCount,Loadi,PD,QD); - %% 开始构建deltF - deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,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,PD,RealPD,RealQD,Loadi,KK,PF); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF,noDataTransCapacity); - 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; +%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; +Volt0=Volt; +UAngel0=UAngel; +%% +PG0(Balance)=PGBal(Balance); +PG(Balance)=PGBal(Balance); +QG0(Balance)=QGBal(Balance); +QG0(PVi)=QGBal(PVi); +QG(PVi)=QGBal(PVi); +%% 真实值 +RealPG=PG0; +RealQG=QG0; +RealPD=PD0; +RealQD=QD0; +%% +[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,RealPD,RealQD,QD,PD); +Gap=(Init_L*Init_Z'-Init_U*Init_W'); +KK=0; +plotGap=zeros(1,60); +ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum; +kmax=60; +Precision=Precision/1; +%% 加误差 +PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1)); +QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1)); +Vbi=sparse(ones(Busnum,1)); +while(abs(Gap)>Precision) + if KK>kmax + break; 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); + 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); %% - %Volt=full(Volt'); - %PD=full(PD); - %% 统计PD误差 - % absPDLoad=abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) ); - absPDLoad=abs( (PD(Loadi)-RealPD(Loadi))./RealPD(Loadi) ); - maxPDError=max(absPDLoad(absPDLoad<10)) - Loadi(maxPDError==absPDLoad) - absQDLoad=abs( (QD(Loadi)-RealQD(Loadi))./RealQD(Loadi) ); - maxQDError=max(absQDLoad(absQDLoad<10)) - Loadi(maxQDError==absQDLoad) - disp('index'); - %Loadi(absPDLoad==maxPDError); - %% 计算总线损 - totalLoss=(sum(PG)-sum(PD(Loadi)))*100; - fprintf('总的损耗为%f(MW 有名值)\n',full(totalLoss)); - fprintf('线损率为 %f\n',full(totalLoss/sum(PG))); - %% 计算各线损 - %Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel); - thesis=thesis.AddPDQDPGQG(PD(Loadi),QD(Loadi),PG(Balance),QG(PVi)); - + 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,RestraintCount,Loadi,PD,QD); + %% 开始构建deltF + deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,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,Vbi); + Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); + Ly=Mat_H; + Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF); + Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,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,Vbi]=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,Vbi); + Gap=(Init_L*Init_Z'-Init_U*Init_W'); + KK=KK+1; end -PD(Loadi)=thesis.MeanPD(); -QD(Loadi)=thesis.MeanQD(); -PG(Balance)=thesis.MeanPG(); -QG(PVi)=thesis.MeanQG(); -thesis.MaxDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi)) -thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi)) -thesis.MaxDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi)) -thesis.StatDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi)) -% thesis.SqareDeviation(RealPG(Balance(1)),RealQG(PVi(1)),RealPD(Loadi),RealQD(Loadi)) -thesis.SqareDeviation(PG0(Balance(1)),QG0(PVi(1)),PD0(Loadi),QD0(Loadi)) -thesis.PercentOfPass(PG0(Balance(1)),QG0(PVi(1)),PD0(Loadi),QD0(Loadi)) -thesis.MaxBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel,'E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle20.txt',PD0,QD0) -thesis.StatBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel) toc 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/OPF_Init.asv b/OPF_Init.asv deleted file mode 100644 index 724c726..0000000 --- a/OPF_Init.asv +++ /dev/null @@ -1,34 +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]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD) -Loadi=find(QD~=0 | PD~=0); -%Loadi=[1:Busnum]'; -RestraintCount=size(Loadi,1)*2+Busnum*1; %约束条件数,放开所有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));% 无功下限 -wPG=0; -wQG=0; -%randInt=randperm(size(Loadi,1)); -%randPDind=randInt(1:10); -randPDind=0; -wPD=1/.05^2*ones(size(Loadi,1),1); -% wPD(1:2:end)=0; -wQD=1/.05^2*zeros(size(Loadi,1),1); -% wQD(1:2:end)=0; -%wD(randPDind)=0;%一些负荷不约束 -%wD(7)=0; -% wD(11)=0; -PD=0.5*PD0; -%powerFacter=0.98; -%QD=PD*sqrt((1-powerFacter^2)/powerFacter^2); -QD=.5*QD0; -end \ No newline at end of file diff --git a/OPF_Init.m b/OPF_Init.m index 6ffb2e4..b475997 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -1,7 +1,7 @@ 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]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD) Loadi=find(QD~=0 | PD~=0); %Loadi=[1:Busnum]'; -RestraintCount=size(Loadi,1)*2+Busnum*1; %约束条件数,放开所有QD +RestraintCount=size(Loadi,1)*2+Busnum*2+Busnum; %约束条件数,放开所有QD t_Bal_volt=Volt(Balance); Volt=sparse(1*ones(1,Busnum)); Volt(Balance)=t_Bal_volt; @@ -20,15 +20,15 @@ wQG=0; %randInt=randperm(size(Loadi,1)); %randPDind=randInt(1:10); randPDind=0; -wPD=1/.05^2*zeros(size(Loadi,1),1); +wPD=zeros(size(Loadi,1),1); wPD(1:2:end)=0; -wQD=1/.05^2*zeros(size(Loadi,1),1); +wQD=zeros(size(Loadi,1),1); wQD(1:2:end)=0; %wD(randPDind)=0;%一些负荷不约束 %wD(7)=0; % wD(11)=0; -PD=0.5*PD0; +PD=0.8*PD0; %powerFacter=0.98; %QD=PD*sqrt((1-powerFacter^2)/powerFacter^2); -QD=.5*QD0; +QD=.8*QD0; end \ No newline at end of file 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/case1V.mat b/case1V.mat deleted file mode 100644 index d303cf9..0000000 Binary files a/case1V.mat and /dev/null differ diff --git a/case2V.mat b/case2V.mat deleted file mode 100644 index e823a8e..0000000 Binary files a/case2V.mat and /dev/null differ diff --git a/case3V.mat b/case3V.mat deleted file mode 100644 index 38606e3..0000000 Binary files a/case3V.mat and /dev/null differ diff --git a/caseM.mat b/caseM.mat deleted file mode 100644 index 1dbbd7f..0000000 Binary files a/caseM.mat and /dev/null differ diff --git a/caseR.mat b/caseR.mat deleted file mode 100644 index 67c6031..0000000 Binary files a/caseR.mat and /dev/null differ diff --git a/func_ddg.asv b/func_ddg.asv deleted file mode 100644 index 43d283b..0000000 --- a/func_ddg.asv +++ /dev/null @@ -1,9 +0,0 @@ -function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD) - -ddg1=sparse(size(PVi,1)+size(PGi,1)+Loadi*2+Busnum,RestraintCount); -t1=sparse(length(Loadi),size(PVi,1)+size(PGi,1)+Busnum) -ddg2=sparse(length(Loadi),diag(-2./PD.^2-2*(QD.^2-PD.^2)/(PD.^2+QD.^2))); -ddg3=diag(2*(PD.^2-QD.^2)/(PD.^2+QD.^2)); -ddg4=sparse(Busnum,RestraintCount); -ddg=[ddg1;ddg2;ddg3;ddg4]; -end \ No newline at end of file diff --git a/func_ddh.m b/func_ddh.m index 7fc87c9..8850a0c 100644 --- a/func_ddh.m +++ b/func_ddh.m @@ -51,6 +51,7 @@ t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT; sizeLoadi=size(Loadi,1)*2; ddh=[ sparse(sizeLoadi,ContrlCount); - sparse(2*Busnum,sizeLoadi),-t; + sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum); + sparse(Busnum,ContrlCount); ]; 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 index a6f63cf..e98aaa6 100644 --- a/func_deltF.m +++ b/func_deltF.m @@ -4,6 +4,7 @@ t4=2*wQD.*(QD(Loadi)-QD0(Loadi)); deltF=[sparse(t3); sparse(t4); sparse(2*Busnum,1); + sparse(Busnum,1); ]; end \ No newline at end of file diff --git a/func_deltG.asv b/func_deltG.asv index 65b99ed..5355449 100644 --- a/func_deltG.asv +++ b/func_deltG.asv @@ -1,46 +1,34 @@ function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD) -%% -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); -%% -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,1); -dg5_dQr=sparse(sizePVi,Busnum); -%% -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),1); +dg4_dPD=sparse(size(Loadi,1),length(Loadi)); dg5_dPD=sparse(size(Loadi,1),Busnum); +dg6_dPD=dg5_dPD; +dg7_dPD=sparse(sizeLoadi,Busnum); %% -dg1_dQD=sparse(size(Loadi(1),1),size(PGi,1)); -dg2_dQD=sparse(size(Loadi(1),1),size(PVi,1)); -dg3_dQD=sparse(length(Loadi(1)),length(Loadi)); -dg4_dQD=sparse(1:size(Loadi(1),1),1:size(Loadi(1),1),ones(size(Loadi(1),1),1),size(Loadi(1),1),size(Loadi(1),1)); -dg5_dQD=sparse(size(Loadi(1),1),Busnum); - +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=dg5_dQD; +dg7_dQD=sparse(sizeLoadi,Busnum); %% -dg1_dx=sparse(2*Busnum,sizePGi); -dg2_dx=sparse(2*Busnum,sizePVi); dg3_dx=sparse(2*Busnum,sizeLoadi); -dg4_dx=sparse(2*Busnum,1); +dg4_dx=sparse(2*Busnum,length(Loadi)); dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum); sparse(Busnum,Busnum); ]; -%% -deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg4_dPg,dg5_dPg; - dg1_dQr,dg2_dQr,dg3_dQr,dg4_dQr,dg5_dQr; - dg1_dPD,dg2_dPD,dg3_dPD,dg4_dPD,dg5_dPD; - dg1_dQD,dg2_dQD,dg3_dQD,dg4_dQD,dg5_dQD; - dg1_dx,dg2_dx,dg3_dx,dg4_dx,dg5_dx; +dg6_dx=dg5_dx; +dg7_dx=sparse(2*Busnum,Busnum); +%% +dg3_dvbi=sparse(Busnum,sizeLoadi); +dg4_dvbi=sparse(Busnum,length(Loadi)); +dg5_dvbi=sparse(Busnum,Busnum); +dg6_dvbi=sparse(Busnum,Busnum); +dg6_dvbi=sparse(Busnum,Busnum); +%% +deltG=[dg3_dPD,dg4_dPD,dg5_dPD,dg6_dPD; + dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD; + dg3_dx,dg4_dx,dg5_dx,dg6_dx; + dg3_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi; ]; end \ No newline at end of file diff --git a/func_deltG.m b/func_deltG.m index 230db27..ee01ea9 100644 --- a/func_deltG.m +++ b/func_deltG.m @@ -1,29 +1,34 @@ function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD) -%% - sizeLoadi=size(Loadi,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=dg5_dPD; +dg7_dPD=sparse(sizeLoadi,Busnum); %% 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=dg5_dQD; +dg7_dQD=sparse(sizeLoadi,Busnum); %% 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); ]; -%% - -deltG=[dg3_dPD,dg4_dPD,dg5_dPD; - dg3_dQD,dg4_dQD,dg5_dQD; - dg3_dx,dg4_dx,dg5_dx; +dg6_dx=dg5_dx; +dg7_dx=sparse(2*Busnum,Busnum); +%% +dg3_dvbi=sparse(Busnum,sizeLoadi); +dg4_dvbi=sparse(Busnum,length(Loadi)); +dg5_dvbi=sparse(Busnum,Busnum); +dg6_dvbi=sparse(Busnum,Busnum); +dg7_dvbi=sparse(eye(Busnum,Busnum)); +%% +deltG=[dg3_dPD,dg4_dPD,dg5_dPD,dg6_dPD,dg7_dPD; + dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD,dg7_dQD; + dg3_dx,dg4_dx,dg5_dx,dg6_dx,dg7_dx; + dg3_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi; ]; 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 index e2b8401..5e411a7 100644 --- a/func_deltH.m +++ b/func_deltH.m @@ -2,6 +2,7 @@ function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi) 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_dPD;dH_dQD;dH_dx']; +dH_dvbi=sparse(Busnum,2*Busnum); +deltH=[dH_dPD;dH_dQD;dH_dx';dH_dvbi]; 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/jacobian.asv b/jacobian.asv deleted file mode 100644 index 1056b37..0000000 --- a/jacobian.asv +++ /dev/null @@ -1,40 +0,0 @@ -function [Jacob,PQ,U,Uangle]=jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c) -%************************************************************************** -% 程序功能 : 子函数——形成雅可比矩阵Jacobian -% 编 者： -% 编制时间：2010.12 -%************************************************************************** -%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q -AngleIJ = Uangle(r) - Uangle(c)- Angle'; -U(PVi) = PVu; -U(Balance)=; -temp1= -sparse(1:Busnum,1:Busnum,U,Busnum,Busnum)*Y*sparse(1:Busnum,1:Busnum,U,Busnum,Busnum); % 计算雅克比矩阵可利用的中间变量 -temp2 = sum(temp1.*sparse(r,c,sin(AngleIJ)),2); -temp3 = sum(temp1.*sparse(r,c,cos(AngleIJ)),2); -temp4=sparse(1:Busnum,1:Busnum,temp2,Busnum,Busnum); -temp5=sparse(1:Busnum,1:Busnum,temp3,Busnum,Busnum); -H = temp1.*sparse(r,c,sin(AngleIJ))-temp4; -L = temp1.*sparse(r,c,sin(AngleIJ))+temp4; -N = temp1.*sparse(r,c,cos(AngleIJ))+temp5; -J = -temp1.*sparse(r,c,cos(AngleIJ))+temp5; - -Q = Q0+temp2'; %求有功分量P -P = P0+temp3'; %求无功分量Q -%% 处理平衡节点和pv节点 -H(:,Balance) = 0; -H(Balance,:) = 0; -%H(Balance,Balance) = 100; % 平衡节点对应的对角元素置一个有限数 -H=H+sparse(Balance,Balance,ones(1,length(Balance)),Busnum,Busnum); -L(:,PVi) = 0; -L(PVi,:) = 0; -L = L+sparse(PVi,PVi,ones(1,length(PVi)),Busnum,Busnum); % PV节点对应的对角元素置为1 -J(:,Balance) = 0; -J(PVi,:) = 0; -N(:,PVi) = 0; -N(Balance,:) = 0; -Q(PVi) = 0; % 将pv节点的无功不平衡分量置零 -P(Balance) = 0; % 平衡节点的有功功率不平衡分量置零 -%% 合成PQ和雅可比矩阵 -PQ = cat(2,P,Q); % 形成功率不平衡分量列向量 -Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % 形成Jacobian矩阵 -end \ No newline at end of file diff --git a/openfile2.asv b/openfile2.asv deleted file mode 100644 index 6d24bb0..0000000 --- a/openfile2.asv +++ /dev/null @@ -1,84 +0,0 @@ -function [Busnum,Balance,PQstandard,Precision,Linei,Linej,Liner,Linex,Lineb,kmax,Transfori ,... - Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL] = openfile2(FileName) -%************************************************************************** -% 程序简介 : 子函数——读取潮流计算所需数据 -% 编 者： -% 编制时间 ：2010.12 -%************************************************************************** -data = dlmread(FileName); % 一次读入全部数据 -zeroRow = find(data(:,1)==0); -Busnum= data(1,1); % 节点数 -PQstandard = data(1,3); % 基准容量 -kmax = data(1,4); %最大迭代次数 -Precision = data(1,4); % 精度 -%Balance = data(3,2); -Balance=data(3:zeroRow(2)-1,2);% 生成1到节点号的列向量 -%CenterA=data(1,5); %中心参数 -%LineNum=data(1,2); %支路数 -Base=data(1,3); -%% 各参数矩阵分块 - -line = data(zeroRow(2)+1:zeroRow(3)-1,:); % 形成线路参数矩阵 -ground = data(zeroRow(5)+1:zeroRow(6)-1,:); % 形成对地支路参数矩阵 -tran = data(zeroRow(3)+1:zeroRow(4)-1,:); % 形成变压器参数矩阵 -buspq = data(zeroRow(8)+1:zeroRow(9)-1,:); % 形成节点功率参数矩阵 -PV = data(zeroRow(11)+1:zeroRow(12)-1,:); % 形成pv节点功率参数矩阵 -GenP=data(zeroRow(9)+1:zeroRow(10)-1,:); -GenQ=data(zeroRow(11)+1:zeroRow(12)-1,:); -%% 线路参数矩阵分块 -Linei = line(:,2); % 节点i -Linej= line(:,3); % 节点j -Liner = line(:,4); % 线路电阻 -Linex = line(:,5); % 线路电抗 -Lineb = line(:,6); % b/2 -%% 对地支路参数矩阵 -Branchi = ground(:,2); % 对地支路节点号 -Branchb = ground(:,4); % 对地支路的导纳 -%% 变压器参数矩阵 -Transfori = tran(:,3); % 节点i -Transforj= tran(:,4); % 节点j -Transforr = tran(:,5); % 变压器电阻 -Transforx= tran(:,6); % 变压器电抗 -Transfork0 = tran(:,7); % 变压器变比 -%% 节点功率参数矩阵 -Pointpoweri = buspq(:,3); -PG=buspq(:,5); % 发电机有功 -QG=buspq(:,6); % 发电机无功 -PD=buspq(:,7); % 负荷有功 -QD=buspq(:,8); % 负荷无功 -%%除以基值 -PG=PG/Base; -QG=QG/Base; -PD=PD/Base; -QD=QD/Base; -%% -PD=sparse(PD); -QD=sparse(QD); -PG=sparse(PG); -QG=sparse(QG); -%% pv节点功率参数矩阵 -PVi = PV(:,3); % PV节点的节点号 -PVu = PV(:,5); % PV节点电压 -PVQL=PV(:,6);%PV节点无功下限 -PVQL=PVQL/Base; -PVQU=PV(:,7); %PV节点无功上限 -PVQU=PVQU/Base; -%% 发电机参数 -%GenU=Gen(:,[1 5 6]); -%GenL=Gen(:,[1 7 8]); -GenC=GenP(:,[3 7:9]); -t=GenC(:,2); -GenC(:,2)=GenC(:,4); -GenC(:,4)=t; -%%%%%%%%%%%%%%%%%%%% -%GenC(:,2:4)=100*GenC(:,2:4); -t=GenP(:,[3 5]); -%GenL=[t,PVQL(PVi)]; -GenL=t;%有功下界 -GenL(:,2)=GenL(:,2)/Base; -t=GenP(:,[3 6]); -%GenU=[t,PVQU(PVi)]; -GenU=t;%有功上届 -GenU(:,2)=GenU(:,2)/Base; -PGi=GenP(:,3);%发电机节点号 -end \ 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/plotVolt.asv b/plotVolt.asv deleted file mode 100644 index 84e68f8..0000000 --- a/plotVolt.asv +++ /dev/null @@ -1,26 +0,0 @@ -%% 画电压，为了写论文用 -load('case1V.mat'); -CV1=Volt; -load('case2V.mat'); -CV2=Volt; -load('case3V.mat'); -CV3=Volt; -load('caseR.mat'); -CVR=Volt; -load('caseM.mat'); -CVM=Volt; -subplot(3,2,1) -hist(CVR) -title('真实值') - - -subplot(2,2,3) -hist(CV1) -title('Case 1') -subplot(2,2,4) -hist(CV2) -title('Case 2') -subplot(2,2,5) -hist(CV3) -title('Case 3') - diff --git a/plotVolt.m b/plotVolt.m index 895db0a..adf5a27 100644 --- a/plotVolt.m +++ b/plotVolt.m @@ -5,6 +5,8 @@ load('case2V.mat'); CV2=Volt; load('case3V.mat'); CV3=Volt; +load('case4V.mat'); +CV4=Volt;%负荷曲线偏差大 load('caseR.mat'); CVR=Volt; load('caseM.mat'); @@ -23,18 +25,22 @@ subplot(3,2,3) hist(CV1) xlabel('电压/p.u'); ylabel('数量/个'); -title('情形1电压分布') +title('100%可知情形电压分布') subplot(3,2,4) hist(CV2) xlabel('电压/p.u'); ylabel('数量/个'); -title('情形2电压分布') +title('50%可知情形电压分布') subplot(3,2,5) hist(CV3) xlabel('电压/p.u'); ylabel('数量/个'); -title('情形3电压分布') - +title('始端覆盖情形电压分布') +subplot(3,2,6) +hist(CV4) +xlabel('电压/p.u'); +ylabel('数量/个'); +title('典型负荷曲线偏差大情形电压分布') % figure() % hold on