diff --git a/OPF.m b/OPF.m index d77ff27..44668c6 100644 --- a/OPF.m +++ b/OPF.m @@ -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'; diff --git a/admmatrix.m b/admmatrix.m index f73ce63..0e462e6 100644 --- a/admmatrix.m +++ b/admmatrix.m @@ -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); \ No newline at end of file diff --git a/func_ddg.m b/func_ddg.m index 9e778fa..bc47135 100644 --- a/func_ddg.m +++ b/func_ddg.m @@ -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 \ No newline at end of file diff --git a/func_ddh3.m b/func_ddh3.m index c710f1a..a9b01ec 100644 --- a/func_ddh3.m +++ b/func_ddh3.m @@ -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 \ No newline at end of file diff --git a/func_deltH.m b/func_deltH.m index 34f7ab8..1c1ea2a 100644 --- a/func_deltH.m +++ b/func_deltH.m @@ -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 需要使用一下转置 暂时改一下 diff --git a/func_deltdeltF.m b/func_deltdeltF.m index 163883d..935d370 100644 --- a/func_deltdeltF.m +++ b/func_deltdeltF.m @@ -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); diff --git a/jacobian_M3.asv b/jacobian_M3.asv index a83f479..87f288b 100644 --- a/jacobian_M3.asv +++ b/jacobian_M3.asv @@ -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 \ No newline at end of file +end \ No newline at end of file diff --git a/jacobian_M3.m b/jacobian_M3.m index 013275a..00053a9 100644 --- a/jacobian_M3.m +++ b/jacobian_M3.m @@ -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 \ No newline at end of file +end \ No newline at end of file diff --git a/openfile.m b/openfile.m index 8171214..7545b7c 100644 --- a/openfile.m +++ b/openfile.m @@ -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节点电压 diff --git a/pf.asv b/pf.asv index 2758c87..575e701 100644 --- a/pf.asv +++ b/pf.asv @@ -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); %求解修正方程,更新电压变量 diff --git a/pf.m b/pf.m index 2758c87..575e701 100644 --- a/pf.m +++ b/pf.m @@ -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); %求解修正方程,更新电压变量