diff --git a/Initial.m b/Initial.m index 71a11df..3ee8ad5 100644 --- a/Initial.m +++ b/Initial.m @@ -8,6 +8,6 @@ function [P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum) P0 = sparse(1, Pointpoweri,(PG-PD)/PQstandard); % 求取节点注入有功功率的标幺值 Q0 = sparse(1, Pointpoweri,(QG-QD)/PQstandard); % 求取节点注入无功功率的标幺值 %% 平启动赋电压初值 -U = ones(1,Busnum); % 按照平启动给电压幅值赋值 +U = 1*ones(1,Busnum); % 按照平启动给电压幅值赋值 Uangle = zeros(1,Busnum); % 按照平启动给电压相角赋值 end \ No newline at end of file diff --git a/OPF.m b/OPF.m index a9a4271..558e1d1 100644 --- a/OPF.m +++ b/OPF.m @@ -2,7 +2,7 @@ tic clc 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,Liner,Linex,Lineb,Transforr,Transforx,Transfork0]= ... -pf('ieee30.dat'); +pf('c:/newFIle.txt'); %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); %pf('c:/file31.txt'); @@ -96,14 +96,15 @@ DrawGap(plotGap); %% 统计PD误差 % absPDLoad=abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) ); absPDLoad=abs( (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi) ); -maxPDError=max(absPDLoad) +maxPDError=max(absPDLoad(absPDLoad<10)) absQDLoad=abs( (QD(Loadi)-QDReal(Loadi))./QDReal(Loadi) ); -maxQDError=max(absQDLoad) +maxQDError=max(absQDLoad(absQDLoad<10)) disp('index'); -Loadi(absPDLoad==maxPDError); +%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); toc diff --git a/admmatrix.m b/admmatrix.m index 0e462e6..16b536f 100644 --- a/admmatrix.m +++ b/admmatrix.m @@ -1,5 +1,5 @@ function [GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori... - ,Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb) + ,Transforj,Transforr,Transforx,Transfork0,Branchi,Branchg,Branchb) %************************************************************************** % 程序功能 : 子函数——形成节点导纳矩阵Y % 编 者： @@ -22,6 +22,7 @@ end %% 接地支路计算 if Branchi>0 % 判断有无接地支路 B = B+sparse(Branchi,Branchi,Branchb,Busnum,Busnum); + G = G+sparse(Branchi,Branchi,Branchg,Busnum,Busnum); end %% 化作极坐标形式 GB = G+B.*1i; %将电导，电纳合并，写成复数形式 diff --git a/jacobian.asv b/jacobian.asv new file mode 100644 index 0000000..1056b37 --- /dev/null +++ b/jacobian.asv @@ -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 \ No newline at end of file diff --git a/jacobian.m b/jacobian.m index c7ecfb9..596b585 100644 --- a/jacobian.m +++ b/jacobian.m @@ -7,6 +7,7 @@ function [Jacob,PQ,U,Uangle]=jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0 %% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q AngleIJ = Uangle(r) - Uangle(c)- Angle'; U(PVi) = PVu; +U(Balance)=1; 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); diff --git a/openfile.m b/openfile.m index 154d8ff..d2afa79 100644 --- a/openfile.m +++ b/openfile.m @@ -51,10 +51,10 @@ QG=QG/Base; PD=PD/Base; QD=QD/Base; %% -PD=sparse(PD); -QD=sparse(QD); -PG=sparse(PG); -QG=sparse(QG); +PD=sparse(PD)/Base; +QD=sparse(QD)/Base; +PG=sparse(PG)/Base; +QG=sparse(QG)/Base; %% pv节点功率参数矩阵 PVi = PV(:,1); % PV节点的节点号 PVu = PV(:,2); % PV节点电压 diff --git a/openfile2.m b/openfile2.m index 69e5250..53952b6 100644 --- a/openfile2.m +++ b/openfile2.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,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL] = openfile2(FileName) + Transforj,Transforr,Transforx,Transfork0,Branchi,Branchg,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL] = openfile2(FileName) %************************************************************************** % 程序简介 : 子函数——读取潮流计算所需数据 % 编 者： @@ -34,6 +34,7 @@ Lineb = line(:,6); % %% 对地支路参数矩阵 Branchi = ground(:,2); % 对地支路节点号 Branchb = ground(:,4); % 对地支路的导纳 +Branchg = ground(:,3); % 对地支路的导纳 %% 变压器参数矩阵 Transfori = tran(:,3); % 节点i Transforj= tran(:,4); % 节点j @@ -54,6 +55,7 @@ QD=QD/Base; %% PD=sparse(PD); QD=sparse(QD); +%QD=PD*sqrt(1-.85^2)/.85; PG=sparse(PG); QG=sparse(QG); %% pv节点功率参数矩阵 diff --git a/pf.m b/pf.m index c3df35a..a1ed312 100644 --- a/pf.m +++ b/pf.m @@ -10,10 +10,10 @@ function [kmax,Precision,Uangle,U,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Li 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); + Transforj,Transforr,Transforx,Transfork0,Branchi,Branchg,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL]= openfile2(FileName); %% 形成节点导纳矩阵 [GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,... - Transforx,Transfork0,Branchi,Branchb); + Transforx,Transfork0,Branchi,Branchg,Branchb); [P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum); %求功率不平衡量 %disp('迭代次数i 最大不平衡量'); %% 循环体计算 @@ -26,13 +26,9 @@ for i = 0:kmax if m > Precision %判断不平衡量是否满足精度要求 [Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U); %求解修正方程，更新电压变量 else - %disp(['收敛，迭代次数为',num2str(i),'次']); + disp(['收敛，迭代次数为',num2str(i),'次']); break %若满足精度要求，则计算收敛 end end toc; -PG=PG/PQstandard; -QG=QG/PQstandard; -PD=PD/PQstandard; -QD=QD/PQstandard; end