facat
/
pes2014-debug2-laplace
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-27 19:54:38 +08:00
parent 97822756be
commit 54cd254cea
12 changed files with 77 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
FormH.m
View File

@ -1,8 +1,8 @@
function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND) function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND)
%% %%
QDcos=textread('1047glys.txt'); % QDcos=textread('1047glys.txt');
QD(QD~=0)=PD(QD~=0)./tan(QDcos); % QD(QD~=0)=PD(QD~=0)./tan(QDcos);
QD(QD_NON_ZERO_IND)=QD_NON_ZERO; % QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
%% %%
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt'; dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';

15
FormLw.asv
View File

@ -1,13 +1,14 @@
function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0) function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK)
PU=1*GenU(:,2);%发电机有功上界 PU=1*GenU(:,2);%发电机有功上界
QU=PVQU(:,1);%发电机无功上界 QU=PVQU(:,1);%发电机无功上界
VoltU=1.2*ones(1,Busnum); VoltU=(1.2+1/ekk)*ones(1,Busnum);
PDU=PD0; PDU=PD0(Loadi);
PDU(PD0>0)=1.2*PD0()'; PDU(PDU>0)=1.2*PDU(PDU>0);
PDU(PDU==0)=20; PDU(PDU<0)=0.8*PDU(PDU<0);
PDU=25*ones(Busnum,1)'; PDU(PDU==0)=.2;
t1=([PU',QU',PDU,VoltU])'; t1=([PU',QU',PDU',VoltU])';
t2=Mat_G+Init_U'-t1; t2=Mat_G+Init_U'-t1;
Lw=t2; Lw=t2;
end end

12
FormLw.m
View File

@ -1,12 +1,12 @@
function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi) function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK)
%KK=999;
PU=1*GenU(:,2);%发电机有功上界 PU=1*GenU(:,2);%发电机有功上界
QU=PVQU(:,1);%发电机无功上界 QU=PVQU(:,1);%发电机无功上界
VoltU=1.1*ones(1,Busnum); VoltU=(1.2+1/exp(KK))*ones(1,Busnum);
PDU=PD0(Loadi); PDU=PD0(Loadi);
PDU(PDU>0)=1.2*PDU(PDU>0); PDU(PDU>0)=1.20*PDU(PDU>0);
PDU(PDU<0)=0.8*PDU(PDU<0); PDU(PDU<0)=.80*PDU(PDU<0);
PDU(PDU==0)=.2; PDU(PDU==0)=.20;
t1=([PU',QU',PDU',VoltU])'; t1=([PU',QU',PDU',VoltU])';
t2=Mat_G+Init_U'-t1; t2=Mat_G+Init_U'-t1;
Lw=t2; Lw=t2;

26
FormLz.asv
View File

@ -1,28 +1,14 @@
function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi) function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK)
%KK=999;
PL=1*GenL(:,2);%发电机有功下界 PL=1*GenL(:,2);%发电机有功下界
QL=PVQL(:,1);%发电机无功下界 QL=PVQL(:,1);%发电机无功下界
VoltL=0.9*ones(1,Busnum); VoltL=(0.8-1/exp(KK))*ones(1,Busnum);
PDL=PD0(Loadi); PDL=PD0(Loadi);
PDL(PD0>0)=0.8*PD0(PDL>0); PDL(PDL>0)=0.8*PDL(PDL>0);
PDL(PD0==0)=-.2; PDL(PDL<0)=1.2*PDL(PDL<0);
PDL(PDL<0)=1.2*PD0(PDL<0); PDL(PDL==0)=-.2;
t1=([PL',QL',PDL',VoltL])'; t1=([PL',QL',PDL',VoltL])';
t2=Mat_G-Init_L'-t1; t2=Mat_G-Init_L'-t1;
Lz=t2; Lz=t2;
% PL=1*GenL(:,2);%发电机有功下界
% QL=PVQL(:,1);%发电机无功下界
% VoltL=0.9*ones(1,Busnum);
% PDL=PD0;
% PDL(PD0>0)=0.8*PD0(PD0>0);
% PDL(PD0==0)=-.2;
% PDL(PD0<0)=1.2*PD0(PD0<0);
%
% t1=([PL',QL',PDL',VoltL])';
%
% t2=Mat_G-Init_L'-t1;
% Lz=t2;
end end

12
FormLz.m
View File

@ -1,12 +1,12 @@
function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi) function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK)
%KK=999;
PL=1*GenL(:,2);%发电机有功下界 PL=1*GenL(:,2);%发电机有功下界
QL=PVQL(:,1);%发电机无功下界 QL=PVQL(:,1);%发电机无功下界
VoltL=0.9*ones(1,Busnum); VoltL=(0.8-1/exp(KK))*ones(1,Busnum);
PDL=PD0(Loadi); PDL=PD0(Loadi);
PDL(PDL>0)=0.8*PDL(PDL>0); PDL(PDL>0)=.80*PDL(PDL>0);
PDL(PDL<0)=1.2*PDL(PDL<0); PDL(PDL<0)=1.20*PDL(PDL<0);
PDL(PDL==0)=-.2; PDL(PDL==0)=-.20;
t1=([PL',QL',PDL',VoltL])'; t1=([PL',QL',PDL',VoltL])';
t2=Mat_G-Init_L'-t1; t2=Mat_G-Init_L'-t1;
Lz=t2; Lz=t2;

22
OPF.asv
View File

@ -1,6 +1,6 @@
tic tic
clear clear
[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL]=pf('ieee10471PG.dat'); [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(c:/iee.txt);
%PVi电压节点序号 %PVi电压节点序号
%PVu电压节点电压标幺值 %PVu电压节点电压标幺值
Volt; Volt;
@ -34,28 +34,28 @@ while(abs(Gap)>Precision)
L_1Z=diag(Init_Z./Init_L); L_1Z=diag(Init_Z./Init_L);
U_1W=diag(Init_W./Init_U); U_1W=diag(Init_W./Init_U);
%% 形成海森阵 %% 形成海森阵
deltdeltF=func_deltdeltF(Busnum,PVi,PGi,wG,wD,ContrlCount); deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount);
%% 形成ddHy %% 形成ddHy
ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi); ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
%% 开始构建ddg %% 开始构建ddg
ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi); ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi);
%% 开始构建deltF %% 开始构建deltF
deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi); deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi);
%% 形成方程矩阵 %% 形成方程矩阵
Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1); Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
Lul=Init_L'.*Init_Z'-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_G=FormG(Volt,PVi,PGi,PG,QG,PD,Loadi);
Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND); Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND);
Ly=Mat_H; Ly=Mat_H;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi); Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi); Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx); 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); 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]=AssignXX1(XX,ContrlCount,RestraintCount,Busnum); [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); [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'); Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=KK+1; KK=KK+1;
@ -64,7 +64,7 @@ fprintf('
ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi) ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
DrawGap(plotGap); DrawGap(plotGap);
%% %%
Volt=full(Volt'); % Volt=full(Volt');
PD=full(PD); % PD=full(PD);
toc toc

12
OPF.m
View File

@ -1,6 +1,6 @@
tic tic
clear clear
[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL]=pf('ieee10471PG.dat'); [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('e:/ËãÀý/standard.txt');
%PVi电压节点序号 %PVi电压节点序号
%PVu电压节点电压标幺值 %PVu电压节点电压标幺值
Volt; Volt;
@ -13,7 +13,7 @@ Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0; KK=0;
plotGap=zeros(1,50); plotGap=zeros(1,50);
ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2; ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
kmax=600; kmax=60;
%% 20120523 临时 %% 20120523 临时
QD_NON_ZERO=QD(PD==0 & QD~=0); QD_NON_ZERO=QD(PD==0 & QD~=0);
QD_NON_ZERO_IND=find(PD==0 & QD~=0); QD_NON_ZERO_IND=find(PD==0 & QD~=0);
@ -48,8 +48,8 @@ while(abs(Gap)>Precision)
Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,Loadi); 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); Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND);
Ly=Mat_H; Ly=Mat_H;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi); Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi); Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
%% 开始解方程 %% 开始解方程
@ -64,7 +64,7 @@ fprintf('
ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi) ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
DrawGap(plotGap); DrawGap(plotGap);
%% %%
Volt=full(Volt'); % Volt=full(Volt');
PD=full(PD); % PD=full(PD);
toc toc

9
OPF_Init.m
View File

@ -3,13 +3,13 @@ Loadi=find(QD~=0);
%Loadi=[1:Busnum]'; %Loadi=[1:Busnum]';
RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1; %约束条件数 RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1; %约束条件数
t_Bal_volt=Volt(Balance); t_Bal_volt=Volt(Balance);
Volt=sparse(ones(1,Busnum)); Volt=sparse(1*ones(1,Busnum));
Volt(Balance)=t_Bal_volt; Volt(Balance)=t_Bal_volt;
UAngel=sparse(1,Busnum); UAngel=sparse(1,Busnum);
Init_Z=sparse(ones(1,RestraintCount)); Init_Z=sparse(ones(1,RestraintCount));
Init_W=sparse(-1*ones(1,RestraintCount)); Init_W=sparse(-1*ones(1,RestraintCount));
Init_L=sparse(ones(1,RestraintCount)); Init_L=1*sparse(ones(1,RestraintCount));
Init_U=sparse(ones(1,RestraintCount)); Init_U=1*sparse(ones(1,RestraintCount));
Init_Y=sparse(1,2*Busnum);%与学姐一致 Init_Y=sparse(1,2*Busnum);%与学姐一致
tPU=sparse(GenU(:,2));% 发电机有功上限 tPU=sparse(GenU(:,2));% 发电机有功上限
tQU=sparse(PVQU(:,1));% 无功上限 tQU=sparse(PVQU(:,1));% 无功上限
@ -19,10 +19,11 @@ PG(PGi)=(tPU+tPL)/2;
QG(PVi)=(tQU+tQL)/2; QG(PVi)=(tQU+tQL)/2;
wG=ones(size(PGi,1),1); wG=ones(size(PGi,1),1);
randInt=randperm(size(Loadi,1)); randInt=randperm(size(Loadi,1));
randPDind=randInt(1:262); randPDind=randInt(1:0);
wD=ones(size(Loadi,1),1); wD=ones(size(Loadi,1),1);
wD(randPDind)=0;%一些负荷不约束 wD(randPDind)=0;%一些负荷不约束
%wD(12)=0;
PD=1*PD0; PD=1*PD0;
end end

26
SolveIt.asv
View File

@ -1,25 +1,31 @@
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) 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); 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'); H=-deltdeltF+ddh;%+ddg*(Init_Z'+Init_W');
t1=diag(Init_L.\Init_Z-Init_U.\Init_W); t1=diag(Init_L.\Init_Z-Init_U.\Init_W);
t2=-deltG*( t1 )*deltG'; t2=-deltG*( t1 )*deltG';
aa=[ aa=[
(H+t2),deltH; (H+t2),deltH;
deltH',zeros(size(Init_Y,),ContrlCount); deltH',zeros(size(Init_Y,2));
]; ];
yy=[LxComa;-Ly]; yy=[LxComa;-Ly];
% t=size(PVi,1)+size(PGi,1); %% ƽºâ½Úµãµçѹ²»±ä
% aa(t+2*Balance-1,:)=0; t=size(PVi,1)+size(PGi,1)+size(Loadi,1);
% aa(:,t+2*Balance-1)=0;
% aa(t+2*Balance-1,t+2*Balance-1)=1;
%ÔÝʱ¸ÄÒ»ÏÂ
t=size(PVi,1)+size(PGi,1)+Busnum;
aa(t+Balance,:)=0; aa(t+Balance,:)=0;
aa(:,t+Balance)=0; aa(:,t+Balance)=0;
aa(t+Balance,t+Balance)=1; %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; dxdy=aa\yy;
dX=dxdy(1:ContrlCount); dX=dxdy(1:ContrlCount);
dY=dxdy(ContrlCount+1:RestraintCount+ContrlCount); dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum);
dL=Lz+deltG'*dX; dL=Lz+deltG'*dX;
dU=-Lw-deltG'*dX; dU=-Lw-deltG'*dX;
dZ=-diag(Init_L)\Lul-diag(Init_L)\diag(Init_Z)*dL; dZ=-diag(Init_L)\Lul-diag(Init_L)\diag(Init_Z)*dL;

6
SolveIt.m
View File

@ -12,13 +12,15 @@ yy=[LxComa;-Ly];
t=size(PVi,1)+size(PGi,1)+size(Loadi,1); t=size(PVi,1)+size(PGi,1)+size(Loadi,1);
aa(t+Balance,:)=0; aa(t+Balance,:)=0;
aa(:,t+Balance)=0; aa(:,t+Balance)=0;
aa(t+Balance,t+Balance)=1; %aa(t+Balance,t+Balance)=1;
aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
deltG(t+Balance,:)=0; deltG(t+Balance,:)=0;
%% %%
t=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1; t=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1;
aa(t+Balance,:)=0; aa(t+Balance,:)=0;
aa(:,t+Balance)=0; aa(:,t+Balance)=0;
aa(t+Balance,t+Balance)=1; %aa(t+Balance,t+Balance)=1;
aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
deltG(t+Balance,:)=0; deltG(t+Balance,:)=0;
%% %%
dxdy=aa\yy; dxdy=aa\yy;

3
jacobian.m
View File

@ -22,7 +22,8 @@ P = P0+temp3';
%% 处理平衡节点和pv节点 %% 处理平衡节点和pv节点
H(:,Balance) = 0; H(:,Balance) = 0;
H(Balance,:) = 0; H(Balance,:) = 0;
H(Balance,Balance) = 100; % 平衡节点对应的对角元素置一个有限数 %H(Balance,Balance) = 100; % 平衡节点对应的对角元素置一个有限数
H=H+sparse(Balance,Balance,ones(1,length(Balance)),Busnum,Busnum);
L(:,PVi) = 0; L(:,PVi) = 0;
L(PVi,:) = 0; L(PVi,:) = 0;
L = L+sparse(PVi,PVi,ones(1,length(PVi)),Busnum,Busnum); % PV节点对应的对角元素置为1 L = L+sparse(PVi,PVi,ones(1,length(PVi)),Busnum,Busnum); % PV节点对应的对角元素置为1

6
openfile.m
View File

@ -6,16 +6,18 @@ function [Busnum,Balance,PQstandard,Precision,Linei,Linej,Liner,Linex,Lineb,kmax
% 编制时间 ：2010.12 % 编制时间 ：2010.12
%************************************************************************** %**************************************************************************
data = dlmread(FileName); % 一次读入全部数据 data = dlmread(FileName); % 一次读入全部数据
zeroRow = find(data(:,1)==0);
Busnum= data(1,1); % 节点数 Busnum= data(1,1); % 节点数
PQstandard = data(1,3); % 基准容量 PQstandard = data(1,3); % 基准容量
kmax = data(1,4); %最大迭代次数 kmax = data(1,4); %最大迭代次数
Precision = data(2,1); % 精度 Precision = data(2,1); % 精度
Balance = data(3,2); % 生成1到节点号的列向量 %Balance = data(3,2);
Balance=data(3:zeroRow(1)-1,2);% 生成1到节点号的列向量
CenterA=data(1,5); %中心参数 CenterA=data(1,5); %中心参数
LineNum=data(1,2); %支路数 LineNum=data(1,2); %支路数
Base=data(1,3); Base=data(1,3);
%% 各参数矩阵分块 %% 各参数矩阵分块
zeroRow = find(data(:,1)==0); %查找第一列元素为零的行号
line = data(zeroRow(1)+1:zeroRow(2)-1,:); % 形成线路参数矩阵 line = data(zeroRow(1)+1:zeroRow(2)-1,:); % 形成线路参数矩阵
ground = data(zeroRow(2)+1:zeroRow(3)-1,:); % 形成对地支路参数矩阵 ground = data(zeroRow(2)+1:zeroRow(3)-1,:); % 形成对地支路参数矩阵
tran = data(zeroRow(3)+1:zeroRow(4)-1,:); % 形成变压器参数矩阵 tran = data(zeroRow(3)+1:zeroRow(4)-1,:); % 形成变压器参数矩阵