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

1.删掉很多没用的文件 

2.基本框架形成，目标函数还没改变。

Signed-off-by: dmy <dugg@21cn.com>
This commit is contained in:
dmy 2014-05-13 16:48:46 +08:00
parent 6fdf435f66
commit 1ee3b9312c
49 changed files with 182 additions and 877 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
0PD0.mat
View File

Binary file not shown.

BIN
0QD0.mat
View File

Binary file not shown.

BIN
12PD0.mat
View File

Binary file not shown.

BIN
12QD0.mat
View File

Binary file not shown.

BIN
16PD0.mat
View File

Binary file not shown.

BIN
16QD0.mat
View File

Binary file not shown.

BIN
20PD0.mat
View File

Binary file not shown.

BIN
20QD0.mat
View File

Binary file not shown.

BIN
4PD0.mat
View File

Binary file not shown.

BIN
4QD0.mat
View File

Binary file not shown.

BIN
8PD0.mat
View File

Binary file not shown.

BIN
8QD0.mat
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;% 太小的值不计算

4
FormG.m
View File

@ -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

8
FormH.m
View File

@ -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';

25
FormLw.asv
View File

@ -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

10
FormLw.m
View File

@ -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;

20
FormLz.asv
View File

@ -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

10
FormLz.m
View File

@ -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;

12
JSMJZM.asv
View File

@ -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')

22
JSMJZM.m
View File

@ -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')

14
JSNZN.asv
View File

@ -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')

38
Lineloss.asv
View File

@ -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

3
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]=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

137
OPF.asv
View File

@ -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

221
OPF.m
View File

@ -4,145 +4,94 @@ 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

106
OPFWrongDataCheck.asv
View File

@ -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

34
OPF_Init.asv
View File

@ -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

10
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)
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

42
SolveIt.asv
View File

@ -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

BIN
case1V.mat
View File

Binary file not shown.

BIN
case2V.mat
View File

Binary file not shown.

BIN
case3V.mat
View File

Binary file not shown.

BIN
caseM.mat
View File

Binary file not shown.

BIN
caseR.mat
View File

Binary file not shown.

9
func_ddg.asv
View File

@ -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

3
func_ddh.m
View File

@ -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

19
func_deltF.asv
View File

@ -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

1
func_deltF.m
View File

@ -4,6 +4,7 @@ t4=2*wQD.*(QD(Loadi)-QD0(Loadi));
deltF=[sparse(t3);
sparse(t4);
sparse(2*Busnum,1);
sparse(Busnum,1);
];
end

54
func_deltG.asv
View File

@ -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);
];
dg6_dx=dg5_dx;
dg7_dx=sparse(2*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;
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

27
func_deltG.m
View File

@ -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);
];
dg6_dx=dg5_dx;
dg7_dx=sparse(2*Busnum,Busnum);
%%
deltG=[dg3_dPD,dg4_dPD,dg5_dPD;
dg3_dQD,dg4_dQD,dg5_dQD;
dg3_dx,dg4_dx,dg5_dx;
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

8
func_deltH.asv
View File

@ -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

3
func_deltH.m
View File

@ -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

11
func_deltdeltF.asv
View File

@ -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

40
jacobian.asv
View File

@ -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

84
openfile2.asv
View File

@ -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

34
pf.asv
View File

@ -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

26
plotVolt.asv
View File

@ -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')

14
plotVolt.m
View File

@ -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