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

优化稀疏化。删掉了很多注释。 

Signed-off-by: facat <dmy@dmy-PC.(none)>
This commit is contained in:
facat 2012-05-24 21:06:34 +08:00
parent 9c5239bd1e
commit f2781393a7
11 changed files with 37 additions and 173 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
OPF.m
View File

@ -1,14 +1,10 @@
tic
clear
%[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,LineLimti,LineLimtj,LinePLimt,PG,QG,PD,QD,CenterA,LineCount,PGi,PVQU,PVQL]=pf('5sj.txt');
[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL]=pf('ieee10471PG.dat');
GB=full(GB);
%PVi电压节点序号
%PVu电压节点电压标幺值
Volt;
UAngel*180/3.1415926;
%sprintf('%f\n',Volt);
%sprintf('%f\n',Angel);
%% 初值
PG0=PG;
PD0=PD;
@ -28,7 +24,6 @@ while(abs(Gap)>Precision)
end
plotGap(KK+1)=Gap;
Init_u=Gap/2/RestraintCount*CenterA;
%AngleIJMat=repmat(UAngel',1,Busnum)-repmat(UAngel,Busnum,1);
AngleIJMat=0;
%% 开始计算OPF
%% 形成等式约束的雅克比
@ -66,7 +61,6 @@ while(abs(Gap)>Precision)
KK=KK+1;
end
fprintf('迭代次数%d\n',KK);
%CalCost(GenC,PG,PGi);
ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD)
DrawGap(plotGap);
Volt=Volt';

4
admmatrix.m
View File

@ -1,4 +1,4 @@
function [G,B,GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori...
function [GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori...
,Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb)
%**************************************************************************
% 程序功能 : 子函数——形成节点导纳矩阵Y
@ -26,7 +26,7 @@ end
%% 化作极坐标形式
GB = G+B.*1i; %将电导，电纳合并，写成复数形式
Y = abs(GB); %求节点导纳幅值
Y=full(Y);
%Y=full(Y);
[r,c] = find(Y);
Angle = angle(GB(GB~=0)); %求节点导纳角度
%Angle=angle(GB);

3
func_ddg.m
View File

@ -1,6 +1,5 @@
function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount)
t=sparse(size(PVi,1)+size(PGi,1)+2*Busnum,RestraintCount);
ddg=sparse(size(PVi,1)+size(PGi,1)+2*Busnum,RestraintCount);
ddg=t;
end

18
func_ddh3.m
View File

@ -7,13 +7,13 @@ mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
mat_INV_AngleIJ=mat_AngleIJ';
yP=Init_Y(1:size(Init_Y,2)/2);%暂时改这里 20111227
yQ=Init_Y(size(Init_Y,2)/2+1:size(Init_Y,2));%暂时改这里 20111227
t1=-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yP');
t2=diag(diag(Volt)*yP')*Y.*cos(mat_AngleIJ);
t3=(t1+t2)*diag(Volt);
t4=-(diag(Y.*cos(mat_AngleIJ)*Volt') -diag(Volt)*Y.*cos(mat_INV_AngleIJ) )*diag(diag(Volt)*yP');
t1=-sparse(1:Busnum,1:Busnum,Y.*cos(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP',Busnum,Busnum);
t2=sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP',Busnum,Busnum)*Y.*cos(mat_AngleIJ);
t3=(t1+t2)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
t4=-(sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt',Busnum,Busnum) -sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP',Busnum,Busnum);
ddPdTdT=t3+t4;%ok1
t1=(-diag(Y.*sin(mat_AngleIJ)*Volt')+diag(Volt)*Y.*sin(mat_INV_AngleIJ) )*diag(yP);
t2= -diag( diag(Volt)*yP' )*Y.*sin(mat_AngleIJ)+diag(Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yP');
t1=(-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt',Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,yP,Busnum,Busnum);
t2= -sparse(1:Busnum,1:Busnum, sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP' ,Busnum,Busnum)*Y.*sin(mat_AngleIJ)+sparse(1:Busnum,1:Busnum,Y.*sin(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP',Busnum,Busnum);
ddPdVdT=t1+t2;%ok1
t1=diag( Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yP');
t2=diag(yP)*Y.*sin(mat_AngleIJ)*diag(Volt);
@ -47,12 +47,8 @@ t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV ;
];
sizePGi=size(PGi,1);
sizePVi=size(PVi,1);
% t=[zeros(size(PGi,1)+size(PVi,1),ContrlCount);
% zeros(2*Busnum,size(PVi,1)+size(PGi,1)),-t;
% ];
t=[
ddh=[
sparse(sizePGi+sizePVi+Busnum,ContrlCount);
sparse(2*Busnum,sizePVi+sizePGi+Busnum),-t;
];
ddh=t;
end

2
func_deltH.m
View File

@ -16,7 +16,7 @@ dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum);
% dH_dQr(I,PVi(I)+Busnum)=1;
% end
dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,1),2*Busnum);
dH_dPD=[-eye(Busnum) zeros(Busnum)];
dH_dPD=[sparse(1:Busnum,1:Busnum,-ones(Busnum,1),Busnum,Busnum) sparse(Busnum,Busnum)];
%Angle=angle(GB);
dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat,UAngel,r,c); %形成雅克比矩阵
%deltH=[dH_dPg;dH_dQr;dH_dx'];%dH_dx 需要使用一下转置 暂时改一下

4
func_deltdeltF.m
View File

@ -4,10 +4,6 @@ ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; %P,Q,Volt theta
C=[wG' zeros(size(PVi))' wD'];
sizeC=size(C,2);
diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC);
% deltdeltF=[diag(C)*2,zeros(size(C,2),ContrlCount-size(C,2));
% zeros(ContrlCount-size(C,2),ContrlCount);
% ];
%sizeGenC=size(GenC(:,2),1);
deltdeltF=[
diagC*2,sparse(sizeC,ContrlCount-sizeC);
sparse(ContrlCount-sizeC,ContrlCount);

113
jacobian_M3.asv
View File

@ -1,4 +1,4 @@
function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat)
function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat,UAngel,r,c)
%**************************************************************************
% 程序功能 : 子函数——形成雅可比矩阵Jacobian
% 编 者：
@ -6,110 +6,15 @@ function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat)
%**************************************************************************
%%参照图书馆6楼的书编写
%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
AngleIJ=AngleIJMat-Angle;
t1=Volt'*Volt;
H=t1.*Y.*sin(AngleIJ);
N=-t1.*Y.*cos(AngleIJ);
%J=Volt'*ones(1,Busnum).*cos(AngleIJ);这里错了
J=Volt'*ones(1,Busnum).*cos(AngleIJ).*Y;
%L=Volt'*ones(1,Busnum).*sin(AngleIJ);这里错了
L=Volt'*ones(1,Busnum).*sin(AngleIJ).*Y;
%%对角
t1=Volt'*Volt;
t2=t1.*Y.*sin(AngleIJ);
t3=diag(t2);
t4=t2-diag(t3);
t5=sum(t4,2);
HH=-diag(t5);
t2=t1.*Y.*cos(AngleIJ);
t3=diag(t2);
t4=t2-diag(t3);
t5=sum(t4,2);
NN=diag(t5);
%t1=ones(Busnum,1)*Volt;
t1=ones(Busnum,1)*Volt.*Y;
t2=t1.*cos(AngleIJ);
t3=sum(t2,2);
JJ=diag(t3);
t1=Volt'*ones(1,Busnum).*cos(AngleIJ).*Y;
%t1=Volt'*ones(1,Busnum).*cos(AngleIJ);
t2=diag(t1);%
JJ=JJ+diag(t2);
t1=ones(Busnum,1)*Volt.*Y;
%t1=ones(Busnum,1)*Volt;
t2=t1.*sin(AngleIJ);
t3=sum(t2,2);
LL=diag(t3);
%t1=Volt'*ones(1,Busnum).*sin(AngleIJ);
t1=Volt'*ones(1,Busnum).*sin(AngleIJ).*Y;
t2=diag(t1);%
%LL=LL-diag(t2);这里错了
LL=LL+diag(t2);
H=H-diag(diag(H));
N=N-diag(diag(N));
J=J-diag(diag(J));
L=L-diag(diag(L));
H=H+HH;
J=J+JJ;
N=N+NN;
L=L+LL;
t1=zeros(2*Busnum);
t1(1:2:2*Busnum,1:2:2*Busnum)=H;
t1(1:2:2*Busnum,2:2:2*Busnum)=N;
t1(2:2:2*Busnum,1:2:2*Busnum)=J;
t1(2:2:2*Busnum,2:2:2*Busnum)=L;
Jacob=-t1;
%%以下是学姐给的公式
H=diag(Volt)*Y.*sin(AngleIJ')*diag(Volt)-diag(Y.*sin(AngleIJ)*Volt')*diag(Volt);
N=-diag(Volt)*Y.*cos(AngleIJ')*diag(Volt)+diag(Y.*cos(AngleIJ)*Volt')*diag(Volt);
J=diag(Y.*cos(AngleIJ)*Volt')+Y.*cos(AngleIJ')*diag(Volt);
L=diag(Y.*sin(AngleIJ)*Volt')+Y.*sin(AngleIJ')*diag(Volt);
H=H;
N=N;
J=J;
L=L;
t1=zeros(2*Busnum);
% t1(1:2:2*Busnum,1:2:2*Busnum)=-H;%10111227
% t1(1:2:2*Busnum,2:2:2*Busnum)=N;
% t1(2:2:2*Busnum,1:2:2*Busnum)=-J;%10111227
% t1(2:2:2*Busnum,2:2:2*Busnum)=L;
%暂时改一下
AngleIJ=UAngel(r)-UAngel(c)-Angle';
mat_AngleIJ=sparse(r,c,AngleIJ,Busnum,Busnum);
mat_IvAngleIJ=mat_AngleIJ';
H=sparse(1:Volt)*Y.*sin(mat_IvAngleIJ)*diag(Volt)-diag(Y.*sin(mat_AngleIJ)*Volt')*diag(Volt);
N=-diag(Volt)*Y.*cos(mat_IvAngleIJ)*diag(Volt)+diag(Y.*cos(mat_AngleIJ)*Volt')*diag(Volt);
J=diag(Y.*cos(mat_AngleIJ)*Volt')+Y.*cos(mat_IvAngleIJ)*diag(Volt);
L=diag(Y.*sin(mat_AngleIJ)*Volt')+Y.*sin(mat_IvAngleIJ)*diag(Volt);
t1=[J,L;
H,N;
]';
Jacob=-t1;
end
% function Jacob=jacobian_M1(Busnum,PVi,PVu,U,Uangle,Y,Angle,r,c)
% AngleIJ = Uangle(r) - Uangle(c)- Angle';
% U(PVi) = PVu;
% 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;
%
%
% t1=zeros(2*Busnum);
% t1(1:2:2*Busnum,1:2:2*Busnum)=H;
% t1(1:2:2*Busnum,2:2:2*Busnum)=N;
% t1(2:2:2*Busnum,1:2:2*Busnum)=J;
% t1(2:2:2*Busnum,2:2:2*Busnum)=L;
% % t1(1:)
% % PQ = cat(2,P,Q); % 形成功率不平衡分量列向量
% %Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % 形成Jacobian矩阵
% Jacob=t1;
%
% end
end

39
jacobian_M3.m
View File

@ -9,41 +9,12 @@ function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat,UAngel,r,c)
AngleIJ=UAngel(r)-UAngel(c)-Angle';
mat_AngleIJ=sparse(r,c,AngleIJ,Busnum,Busnum);
mat_IvAngleIJ=mat_AngleIJ';
H=diag(Volt)*Y.*sin(mat_IvAngleIJ)*diag(Volt)-diag(Y.*sin(mat_AngleIJ)*Volt')*diag(Volt);
N=-diag(Volt)*Y.*cos(mat_IvAngleIJ)*diag(Volt)+diag(Y.*cos(mat_AngleIJ)*Volt')*diag(Volt);
J=diag(Y.*cos(mat_AngleIJ)*Volt')+Y.*cos(mat_IvAngleIJ)*diag(Volt);
L=diag(Y.*sin(mat_AngleIJ)*Volt')+Y.*sin(mat_IvAngleIJ)*diag(Volt);
H=sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt',Busnum,Busnum)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
N=-sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt',Busnum,Busnum)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
J=sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt',Busnum,Busnum)+Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
L=sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt',Busnum,Busnum)+Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
t1=[J,L;
H,N;
]';
Jacob=-t1;
end
% function Jacob=jacobian_M1(Busnum,PVi,PVu,U,Uangle,Y,Angle,r,c)
% AngleIJ = Uangle(r) - Uangle(c)- Angle';
% U(PVi) = PVu;
% 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;
%
%
% t1=zeros(2*Busnum);
% t1(1:2:2*Busnum,1:2:2*Busnum)=H;
% t1(1:2:2*Busnum,2:2:2*Busnum)=N;
% t1(2:2:2*Busnum,1:2:2*Busnum)=J;
% t1(2:2:2*Busnum,2:2:2*Busnum)=L;
% % t1(1:)
% % PQ = cat(2,P,Q); % 形成功率不平衡分量列向量
% %Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % 形成Jacobian矩阵
% Jacob=t1;
%
% end
end

7
openfile.m
View File

@ -1,5 +1,5 @@
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)
Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL] = openfile(FileName)
%**************************************************************************
% 程序简介 : 子函数——读取潮流计算所需数据
% 编 者：
@ -48,6 +48,11 @@ 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(:,1); % PV节点的节点号
PVu = PV(:,2); % PV节点电压

7
pf.asv
View File

@ -10,9 +10,9 @@ 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,Gen,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL]= openfile(FileName);
Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL]= openfile(FileName);
%% 形成节点导纳矩阵
[G,B,GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,...
[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 最大不平衡量');
@ -20,9 +20,8 @@ disp('
for i = 0:kmax
[Jacob,PQ,U,Uangle] = jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c); %形成雅克比矩阵
% disp('第一次雅克比');
%full(Jacob);
m = max(abs(PQ));
m = full(m);
m=full(m);
fprintf(' %u %.8f \n',i,m);
if m > Precision %判断不平衡量是否满足精度要求
[Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U); %求解修正方程，更新电压变量

7
pf.m
View File

@ -10,9 +10,9 @@ 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,Gen,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL]= openfile(FileName);
Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL]= openfile(FileName);
%% 形成节点导纳矩阵
[G,B,GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,...
[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 最大不平衡量');
@ -20,9 +20,8 @@ disp('
for i = 0:kmax
[Jacob,PQ,U,Uangle] = jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c); %形成雅克比矩阵
% disp('第一次雅克比');
%full(Jacob);
m = max(abs(PQ));
m = full(m);
m=full(m);
fprintf(' %u %.8f \n',i,m);
if m > Precision %判断不平衡量是否满足精度要求
[Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U); %求解修正方程，更新电压变量