facat/pes2014
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Compare commits

base: facat:20140513
Branches Tags
facat:master facat:20140513 facat:thesis facat:withVolt facat:dev
..
compare: facat:withVolt
Branches Tags
facat:master facat:20140513 facat:thesis facat:withVolt facat:dev

1 Commits
20140513 ... withVolt

Author SHA1 Message Date
facat
d03434dcc5 将电压加入目标函数中。 2013-04-04 22:50:20 +08:00
61 changed files with 1131 additions and 250 deletions
Expand all files Collapse all files

1
.gitignore vendored
View File

@@ -1 +0,0 @@
*.asv

BIN
0PD0.mat Normal file
View File

Binary file not shown.

BIN
0QD0.mat Normal file
View File

Binary file not shown.

BIN
12PD0.mat Normal file
View File

Binary file not shown.

BIN
12QD0.mat Normal file
View File

Binary file not shown.

BIN
16PD0.mat Normal file
View File

Binary file not shown.

BIN
16QD0.mat Normal file
View File

Binary file not shown.

BIN
20PD0.mat Normal file
View File

Binary file not shown.

BIN
20QD0.mat Normal file
View File

Binary file not shown.

BIN
4PD0.mat Normal file
View File

Binary file not shown.

BIN
4QD0.mat Normal file
View File

Binary file not shown.

BIN
8PD0.mat Normal file
View File

Binary file not shown.

BIN
8QD0.mat Normal file
View File

Binary file not shown.

4
@ForThesis/MaxBranchDeviation.m
View File

@@ -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;% 太小的值不计算

6
DrawGap.asv Normal file
View File

@@ -0,0 +1,6 @@
function DrawGap(plotGap)
x=find(plotGap);
ts=size(x,2);
end

12
FormAA.asv Normal file
View File

@@ -0,0 +1,12 @@
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

7
FormAlphaP.asv Normal file
View File

@@ -0,0 +1,7 @@
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

11
FormG.asv Normal file
View File

@@ -0,0 +1,11 @@
function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi)
Mat_G=[
PG(PGi);
QG(PVi);
sparse(PD(Loadi));
sparse(QD(Loadi));
Volt';
];
end

14
FormG.m
View File

@@ -1,14 +1,8 @@
function Mat_G=FormG(Volt,PD,QD,Loadi,Vbi,mVolt,RealVolt,bigM,sigma,PDbi,QDbi,RealPD,RealQD)
function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi)
Mat_G=[
sparse(PD(Loadi))-bigM*PDbi-RealPD(Loadi)*sigma;
sparse(PD(Loadi))+bigM*PDbi+RealPD(Loadi)*sigma;
sparse(QD(Loadi))-bigM*QDbi-RealQD(Loadi)*sigma;
sparse(QD(Loadi))+bigM*QDbi+RealQD(Loadi)*sigma;
Volt'-bigM*Vbi-mVolt'-sigma*RealVolt';
Volt'+bigM*Vbi+mVolt'+sigma*RealVolt';
Vbi;
PDbi;
QDbi;
sparse(PD(Loadi));
sparse(QD(Loadi));
Volt';
];
end

31
FormH.asv Normal file
View File

@@ -0,0 +1,31 @@
function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y)
t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=Volt'*Volt;
t3=t1.*t2;
t4=sum(-t3,2);%P
t5=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t6=t2.*t5;
t7=sum(-t6,2);%Q
t8=PG-PD;
t9=QG-QD;
%Mat_H=([(PG-PD)',(QG-QD)'])'+([t4',t7'])';
Mat_H(1:2:2*Busnum)=t8(1:Busnum)+t4(1:Busnum);
Mat_H(2:2:2*Busnum)=t9(1:Busnum)+t7(1:Busnum);
Mat_H=Mat_H';
%%
QDcos=textread('300glys.txt');
t=QD(PD==0 &&);
aa=QD;
QD(QD~=0)=PD(QD~=0)./tan(QDcos);
QD(PD==0)=t;
%%
%%%%一下是学姐给的公式
AngleIJ=AngleIJMat-angle(GB);
%dP=PG-PD-diag(Volt)*Y*cos(AngleIJ)*Volt';
dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';
%Mat_H(1:2:2*Busnum)=dP(1:Busnum);暂时改一下 20111227
%Mat_H(2:2:2*Busnum)=dQ(1:Busnum);暂时改一下 20111227
Mat_H=[dP;dQ;];
end

8
FormH.m
View File

@@ -1,4 +1,10 @@
function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi)
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));
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';

25
FormLw.asv Normal file
View File

@@ -0,0 +1,25 @@
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

20
FormLw.m
View File

@@ -1,24 +1,22 @@
function [Lw,UpperLimit]=FormLw(Mat_G,Init_U,Busnum,PD0,QD0,Loadi,bigM)
function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity)
KK=999;
%PU=GenU(:,2);%·¢µç»úÓй¦ÉϽç
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;
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));
t1=([0*PDU',bigM*ones(1,length(Loadi)),0*QDU',bigM*ones(1,length(Loadi)),0*VoltU,bigM*ones(1,Busnum),1*ones(1,Busnum),1*ones(1,length(Loadi)*2)])';
UpperLimit=t1;
% PF=0.85;
% QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
t1=([PDU',QDU',VoltU])';
t2=Mat_G+Init_U'-t1;
Lw=t2;

20
FormLz.asv Normal file
View File

@@ -0,0 +1,20 @@
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

20
FormLz.m
View File

@@ -1,24 +1,20 @@
function [Lz,LoweLimit]=FormLz(Mat_G,Init_L,Busnum,PD0,QD0,Loadi,bigM)
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.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));
%t1=([-1*ones(1,length(Loadi)),0*PDL',-1*ones(1,length(Loadi)),0*QDL',-5*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum),0*ones(1,length(Loadi)*2)])';
t1=([-bigM*ones(1,length(Loadi)),0*PDL',-bigM*ones(1,length(Loadi)),0*QDL',-bigM*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum),0*ones(1,length(Loadi)*2)])';
LoweLimit=t1;
% QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
t1=([PDL',QDL',VoltL])';
t2=Mat_G-Init_L'-t1;
Lz=t2;
end

5
FormYY.asv Normal file
View File

@@ -0,0 +1,5 @@
function FormYY(Init_L,Lul.Lz,Init_U,Luu,Lw,Lz,LxComa)
t=[
-
]
end

12
JSMJZM.asv Normal file
View File

@@ -0,0 +1,12 @@
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')

22
JSMJZM.m
View File

@@ -1,18 +1,18 @@
JSM=[0.070244 0.0533 0.004379];
JZM=[0.002275 0.047311 0.070318];
for I=0.01:0.01:0.08
JSM=[0.0575 0.051 0.0029];
JZM=[0.0003 0.02629 0.0582];
for I=0.01:0.01:0.06
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.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')
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')
hold on
plot(1:3,JZM,'k--');
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')
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')

14
JSNZN.asv Normal file
View File

@@ -0,0 +1,14 @@
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')

38
Lineloss.asv Normal file
View File

@@ -0,0 +1,38 @@
%% 计算线损
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

18
Modification.asv Normal file
View File

@@ -0,0 +1,18 @@
function [Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi)
AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
fprintf('AlphaP %f\n',AlphaP);
AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
fprintf('AlphaD %f\n',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(PVi)=PG(PGi)+deltX(size(PGi,1));
QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1,1:(size(PVi,1)+size(PGi,1)));
t=deltX(size(PVi,1)+size(PGi,1)+1,ContrlCount)';
t(2*Balance-1)=0;
Volt=Volt+AlphaP*t(2:2:2*Busnum);
UAngel=UAngel+AlphaP*t(1:2:2*Busnum);
end

5
Modification.m
View File

@@ -1,4 +1,4 @@
function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=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,PDbi,QDbi)
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);
@@ -24,7 +24,4 @@ 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:2*Busnum+Busnum)';
PDbi=PDbi+AlphaP*t(2*Busnum+Busnum+1:2*Busnum+Busnum+length(Loadi))';
QDbi=QDbi+AlphaP*t(2*Busnum+Busnum+length(Loadi)+1:end)';
end

137
OPF.asv Normal file
View File

@@ -0,0 +1,137 @@
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

236
OPF.m
View File

@@ -3,108 +3,146 @@ clc
clear
%% 存在问题
% 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
%%
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('E:\算例\东际911_2751267_2012-09-05\newFIle9-2.txt');
%%
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');
% 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;
%%
[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+length(Loadi)*2;
kmax=600;
Precision=Precision/1;
sigma=0.03;
%% 加误差
PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
mPD=PD0;
mQD=QD0;
RealVolt=Volt0;
mVolt=Volt0.*(1+normrnd(0,sigma,1,length(Volt0)));
Vbi=sparse(1*ones(Busnum,1));
PDbi=sparse(1*ones(length(Loadi),1));
QDbi=sparse(1*ones(length(Loadi),1));
bigM=5;
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);
%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))
%%
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,PD,QD,Loadi,Vbi,mVolt,RealVolt,bigM,sigma,PDbi,QDbi,RealPD,RealQD);
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
Ly=Mat_H;
[Lz,LowerLimit]=FormLz(Mat_G,Init_L,Busnum,PD0,QD0,Loadi,bigM);%下界
[Lw,UpperLimit]=FormLw(Mat_G,Init_U,Busnum,PD0,QD0,Loadi,bigM);%上界
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
if any(LowerLimit>UpperLimit)
warning('不等式设置错误');
end
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,PDbi,QDbi]=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,PDbi,QDbi);
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,wVolt,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=KK+1;
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,wVolt,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(Volt,PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,wVolt,PG0,QG0,PD0,PD,QD,QD0,Volt0,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
fprintf('目标函数: %f\n',sum(full(Vbi)));
fprintf('目标函数: %f\n',sum(full(QDbi)));
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

106
OPFWrongDataCheck.asv Normal file
View File

@@ -0,0 +1,106 @@
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

34
OPF_Init.asv Normal file
View File

@@ -0,0 +1,34 @@
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

17
OPF_Init.m
View File

@@ -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)
function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,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+size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2; %约束条件数,放开所有QD
RestraintCount=size(Loadi,1)*2+Busnum*1; %约束条件数,放开所有QD
t_Bal_volt=Volt(Balance);
Volt=sparse(1*ones(1,Busnum));
Volt(Balance)=t_Bal_volt;
@@ -20,15 +20,16 @@ wQG=0;
%randInt=randperm(size(Loadi,1));
%randPDind=randInt(1:10);
randPDind=0;
wPD=zeros(size(Loadi,1),1);
wPD(1:2:end)=0;
wQD=zeros(size(Loadi,1),1);
wQD(1:2:end)=0;
wPD=1/.05^2*ones(size(Loadi,1),1);
% wPD(1:2:end)=0;
wQD=1/.05^2*ones(size(Loadi,1),1);
% wQD(1:2:end)=0;
%wD(randPDind)=0;%一些负荷不约束
%wD(7)=0;
% wD(11)=0;
PD=0.8*PD0;
PD=0.5*PD0;
%powerFacter=0.98;
%QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
QD=.8*QD0;
QD=.5*QD0;
wVolt=1/.05^2*ones(size(Busnum,1),1);
end

42
SolveIt.asv Normal file
View File

@@ -0,0 +1,42 @@
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

31
admmatrix.asv Normal file
View File

@@ -0,0 +1,31 @@
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);

BIN
case1V.mat Normal file
View File

Binary file not shown.

BIN
case2V.mat Normal file
View File

Binary file not shown.

BIN
case3V.mat Normal file
View File

Binary file not shown.

BIN
caseM.mat Normal file
View File

Binary file not shown.

BIN
caseR.mat Normal file
View File

Binary file not shown.

9
func_ddg.asv Normal file
View File

@@ -0,0 +1,9 @@
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

3
func_ddh.m
View File

@@ -51,7 +51,6 @@ t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
sizeLoadi=size(Loadi,1)*2;
ddh=[
sparse(sizeLoadi,ContrlCount);
sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum+length(Loadi)*2);
sparse(Busnum+length(Loadi)*2,ContrlCount);
sparse(2*Busnum,sizeLoadi),-t;
];
end

19
func_deltF.asv Normal file
View File

@@ -0,0 +1,19 @@
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

14
func_deltF.m
View File

@@ -1,11 +1,11 @@
function deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi)
function deltF=func_deltF(Volt,PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,wVolt,PG0,QG0,PD0,PD,QD,QD0,Volt0,Busnum,Loadi)
t3=2*wPD.*(PD(Loadi)-PD0(Loadi));
t4=2*wQD.*(QD(Loadi)-QD0(Loadi));
deltF=[sparse(length(Loadi),1);
sparse(length(Loadi),1);
sparse(2*Busnum,1);
sparse(ones(Busnum,1));
sparse(ones(length(Loadi),1));
sparse(ones(length(Loadi),1));
t5=2*wVolt.*(Volt-Volt0);
deltF=[sparse(t3);
sparse(t4);
sparse(t5');
sparse(1*Busnum,1);
];
end

46
func_deltG.asv Normal file
View File

@@ -0,0 +1,46 @@
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);
dg5_dPD=sparse(size(Loadi,1),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);
%%
dg1_dx=sparse(2*Busnum,sizePGi);
dg2_dx=sparse(2*Busnum,sizePVi);
dg3_dx=sparse(2*Busnum,sizeLoadi);
dg4_dx=sparse(2*Busnum,1);
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;
];
end

78
func_deltG.m
View File

@@ -1,83 +1,29 @@
function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
%%
sizeLoadi=size(Loadi,1);
%dg3 PD-M*b-t-PD0
%%
%%
dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
%dg32 PD+M*B+t+PD0
dg32_dPD=dg3_dPD;
%dg4 QD-M*B-t-QD0
dg4_dPD=sparse(size(Loadi,1),length(Loadi));
%dg42 QD+M*B+t+QD0
dg42_dPD=sparse(size(Loadi,1),length(Loadi));
dg5_dPD=sparse(size(Loadi,1),Busnum);
dg6_dPD=dg5_dPD;
%dg7 ÊÇ VbiµÄÔ¼Êø
dg7_dPD=sparse(sizeLoadi,Busnum);
dpdbi_dPD=sparse(sizeLoadi,sizeLoadi);
dqdbi_dPD=sparse(sizeLoadi,sizeLoadi);
%%
dg3_dQD=sparse(length(Loadi),length(Loadi));
dg32_dQD=sparse(sizeLoadi,sizeLoadi);
dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
dg42_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);
dpdbi_dQD=sparse(sizeLoadi,sizeLoadi);
dqdbi_dQD=sparse(sizeLoadi,sizeLoadi);
%%
dg3_dx=sparse(2*Busnum,sizeLoadi);
dg32_dx=sparse(2*Busnum,sizeLoadi);
dg4_dx=sparse(2*Busnum,length(Loadi));
dg42_dx=sparse(2*Busnum,sizeLoadi);
dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
sparse(Busnum,Busnum);
];
dg6_dx=dg5_dx;
dg7_dx=sparse(2*Busnum,Busnum);
dpdbi_dx=sparse(2*Busnum,sizeLoadi);
dqdbi_dx=sparse(2*Busnum,sizeLoadi);
%%
dg3_dvbi=sparse(Busnum,sizeLoadi);
dg32_dvbi=sparse(Busnum,sizeLoadi);
dg4_dvbi=sparse(Busnum,length(Loadi));
dg42_dvbi=sparse(Busnum,sizeLoadi);
dg5_dvbi=sparse(-eye(Busnum,Busnum));
dg6_dvbi=sparse(eye(Busnum,Busnum));
dg7_dvbi=sparse(eye(Busnum,Busnum));
dpdbi_dvbi=sparse(Busnum,sizeLoadi);
dqdbi_dvbi=sparse(Busnum,sizeLoadi);
%%
dg3_dpdbi=sparse(-eye(sizeLoadi,sizeLoadi));
dg32_dpdbi=sparse(eye(sizeLoadi,sizeLoadi));
dg4_dpdbi=sparse(sizeLoadi,sizeLoadi);
dg42_dpdbi=sparse(sizeLoadi,sizeLoadi);
dg5_dpdbi=sparse(sizeLoadi,Busnum);
dg6_dpdbi=sparse(sizeLoadi,Busnum);
dg7_dpdbi=sparse(sizeLoadi,Busnum);
dpdbi_dpdbi=sparse(eye(sizeLoadi,sizeLoadi));
dqdbi_dpdbi=sparse(sizeLoadi,sizeLoadi);
%%
dg3_dqdbi=sparse(sizeLoadi,sizeLoadi);
dg32_dqdbi=sparse(sizeLoadi,sizeLoadi);
dg4_dqdbi=sparse(-eye(sizeLoadi,sizeLoadi));
dg42_dqdbi=sparse(eye(sizeLoadi,sizeLoadi));
dg5_dqdbi=sparse(sizeLoadi,Busnum);
dg6_dqdbi=sparse(sizeLoadi,Busnum);
dg7_dqdbi=sparse(sizeLoadi,Busnum);
dpdbi_dqdbi=sparse(sizeLoadi,sizeLoadi);
dqdbi_dqdbi=sparse(eye(sizeLoadi,sizeLoadi));
%%
deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg42_dPD,dg5_dPD,dg6_dPD,dg7_dPD,dpdbi_dPD,dqdbi_dPD;
dg3_dQD,dg32_dQD,dg4_dQD,dg42_dQD,dg5_dQD,dg6_dQD,dg7_dQD,dpdbi_dQD,dqdbi_dQD;
dg3_dx,dg32_dx,dg4_dx,dg42_dx,dg5_dx,dg6_dx,dg7_dx,dpdbi_dx,dqdbi_dx;
dg3_dvbi,dg32_dvbi,dg4_dvbi,dg42_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi,dpdbi_dvbi,dqdbi_dvbi;
dg3_dpdbi,dg32_dpdbi,dg4_dpdbi,dg42_dpdbi,dg5_dpdbi,dg6_dpdbi,dg7_dpdbi,dpdbi_dpdbi,dqdbi_dpdbi;
dg3_dqdbi,dg32_dqdbi,dg4_dqdbi,dg42_dqdbi,dg5_dqdbi,dg6_dqdbi,dg7_dqdbi,dpdbi_dqdbi,dqdbi_dqdbi;
%%
deltG=[dg3_dPD,dg4_dPD,dg5_dPD;
dg3_dQD,dg4_dQD,dg5_dQD;
dg3_dx,dg4_dx,dg5_dx;
];
% deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg5_dPD,dg6_dPD,dg7_dPD;
% dg3_dQD,dg32_dQD,dg4_dQD,dg5_dQD,dg6_dQD,dg7_dQD;
% dg3_dx,dg32_dx,dg4_dx,dg5_dx,dg6_dx,dg7_dx;
% dg3_dvbi,dg32_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi;
% dg3_dpdbi,dg32_dpdbi,dg4_dpdbi,dg5_dpdbi,dg6_dpdbi,dg7_dpdbi;
% ];
end

8
func_deltH.asv Normal file
View File

@@ -0,0 +1,8 @@
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

5
func_deltH.m
View File

@@ -2,9 +2,6 @@ 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); %形成雅克比矩阵
dH_dvbi=sparse(Busnum,2*Busnum);
dH_dpdbi=sparse(length(Loadi),2*Busnum);
dH_dqdbi=sparse(length(Loadi),2*Busnum);
deltH=[dH_dPD;dH_dQD;dH_dx';dH_dvbi;dH_dpdbi;dH_dqdbi];
deltH=[dH_dPD;dH_dQD;dH_dx'];
end

11
func_deltdeltF.asv Normal file
View File

@@ -0,0 +1,11 @@
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

12
func_deltdeltF.m
View File

@@ -1,3 +1,11 @@
function deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount)
deltdeltF=0;
function deltdeltF=func_deltdeltF(PVi,wPG,wQG,wVolt,wPD,wQD,ContrlCount)
%ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数
C=[wPD' wQD' wVolt];
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

40
jacobian.asv Normal file
View File

@@ -0,0 +1,40 @@
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

18
modifyadmmatrix.asv Normal file
View File

@@ -0,0 +1,18 @@
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)); %求节点导纳角度

64
openfile.asv Normal file
View File

@@ -0,0 +1,64 @@
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,Gen,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL] = openfile(FileName)
%**************************************************************************
% 程序简介 : 子函数——读取潮流计算所需数据
% 编 者：
% 编制时间 ：2010.12
%**************************************************************************
data = dlmread(FileName); % 一次读入全部数据
Busnum= data(1,1); % 节点数
PQstandard = data(1,3); % 基准容量
kmax = data(1,4); %最大迭代次数
Precision = data(2,1); % 精度
Balance = data(3,2); % 生成1到节点号的列向量
CenterA=data(1,5); %中心参数
LineNum=data(1,2); %支路数
%% 各参数矩阵分块
zeroRow = find(data(:,1)==0); %查找第一列元素为零的行号
line = data(zeroRow(1)+1:zeroRow(2)-1,:); % 形成线路参数矩阵
ground = data(zeroRow(2)+1:zeroRow(3)-1,:); % 形成对地支路参数矩阵
tran = data(zeroRow(3)+1:zeroRow(4)-1,:); % 形成变压器参数矩阵
buspq = data(zeroRow(4)+1:zeroRow(5)-1,:); % 形成节点功率参数矩阵
PV = data(zeroRow(5)+1:zeroRow(6)-1,:); % 形成pv节点功率参数矩阵
Gen=data(zeroRow(6)+1:zeroRow(7)-1,:);
%% 线路参数矩阵分块
Linei = line(:,2); % 节点i
Linej= line(:,3); % 节点j
Liner = line(:,4); % 线路电阻
Linex = line(:,5); % 线路电抗
Lineb = line(:,6); % b/2
%% 对地支路参数矩阵
Branchi = ground(:,1); % 对地支路节点号
Branchb = ground(:,2); % 对地支路的导纳
%% 变压器参数矩阵
Transfori = tran(:,2); % 节点i
Transforj= tran(:,3); % 节点j
Transforr = tran(:,4); % 变压器电阻
Transforx= tran(:,5); % 变压器电抗
Transfork0 = tran(:,6); % 变压器变比
%% 节点功率参数矩阵
Pointpoweri = buspq(:,1);
PG=buspq(:,2); % 发电机有功
QG=buspq(:,3); % 发电机无功
PD=buspq(:,4); % 负荷有功
QD=buspq(:,5); % 负荷无功
%% pv节点功率参数矩阵
PVi = PV(:,1); % PV节点的节点号
PVu = PV(:,2); % PV节点电压
PVQL=PV(:,3);%PV节点无功下限
PVQU=PV(:,4); %PV节点无功上限
%% 发电机参数
%GenU=Gen(:,[1 5 6]);
%GenL=Gen(:,[1 7 8]);
GenC=Gen(:,[1 2:4]);
t=GenC(:,2);
GenC(:,2)=GenC(:,4);
GenC(:,4)=t;
t=Gen(:,[1 5]);
%GenL=[t,PVQL(PVi)];
GenL=t;%有功下界
t=Gen(:,[1 6]);
%GenU=[t,PVQU(PVi)];
GenU=t;%有功上届
PGi=Gen(:,1);%发电机节点号
end

84
openfile2.asv Normal file
View File

@@ -0,0 +1,84 @@
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

34
pf.asv Normal file
View File

@@ -0,0 +1,34 @@
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

26
plotVolt.asv Normal file
View File

@@ -0,0 +1,26 @@
%% 画电压，为了写论文用
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')

14
plotVolt.m
View File

@@ -5,8 +5,6 @@ load('case2V.mat');
CV2=Volt;
load('case3V.mat');
CV3=Volt;
load('case4V.mat');
CV4=Volt;%负荷曲线偏差大
load('caseR.mat');
CVR=Volt;
load('caseM.mat');
@@ -25,22 +23,18 @@ subplot(3,2,3)
hist(CV1)
xlabel('电压/p.u');
ylabel('数量/个');
title('100%可知情形电压分布')
title('情形1电压分布')
subplot(3,2,4)
hist(CV2)
xlabel('电压/p.u');
ylabel('数量/个');
title('50%可知情形电压分布')
title('情形2电压分布')
subplot(3,2,5)
hist(CV3)
xlabel('电压/p.u');
ylabel('数量/个');
title('始端覆盖情形电压分布')
subplot(3,2,6)
hist(CV4)
xlabel('电压/p.u');
ylabel('数量/个');
title('典型负荷曲线偏差大情形电压分布')
title('情形3电压分布')
% figure()
% hold on