%******************************************************* % 程序名称：电力系统无功优化程序 % 程序功能：主程序 % 程序算法：现代内点法，采用NCP函数处理离散变量 % 目标函数：发电费用最小（1）、有功出力最小（2）、有功网损最小（3）、无功网损最小（4） % 编者： 韦园清 % 编写时间：2010年11月 %******************************************************* function []=OPF(filename) clc; %% 读取数据 data=dlmread(filename); if data(1,1)<1780 %读取1780节点系统格式类型的数据 [nodeNum,lineNum,capacity,kmax,centrPara,accuracy,opfGoal,balNum,balNode,mtrLine,mtrGround,mtrTrans,runPara,Pg,Qg,Pd,Qd,pvNum,pvNode,pvV,... pvQmin,pvQmax,pgNumall,pgNodeall,pgNum,pgNode,a,b,c,Pmin,Pmax,capNum,capI,capGi,capBi,capK,capKmin,capKmax,Vmin,Vmax]=readData(data); tic; %形成节点导纳矩阵 [G,B,transY,transG,transB,transG1,transB1,transNum,lineNum,lineI,lineJ,lineB,transI,transJ,transK0,transKmin,transKmax]=formY(nodeNum,... mtrLine,mtrGround,mtrTrans,capI,capK); else %读取IEEE标准格式数据 [nodeNum,lineNum,capacity,kmax,accuracy,balNum,balNode,mtrLine,mtrGround,mtrTrans,runPara,Pg,Qg,Pd,Qd,pvNum,pvNode,pvV,pvQmin,pvQmax,... mtrCapacitor,pgNum,pgNode,a,b,c,Pmin,Pmax,Vmin,Vmax,centrPara]=readData_RPP(data) ; tic; %形成节点导纳矩阵 [G,B,lineI,lineJ,lineB,transI,transJ,transK0,transNum,transKmin,transKmax,capI,capGi,capBi,capK,capKmin,capKmax,capNum,transG,transB,... transG1,transB1,lineNum]=formY_RPP(nodeNum,lineNum,mtrLine,mtrGround,mtrTrans,mtrCapacitor) ; end %% 初始化 [e,f,PG,QR,Pg,Qg,T,pgvNum,xNum,m,r,l,u,z,w,x,y,mu,dmaxOut1,dmaxOut2,xz1,xz2,dmaxOut3,dmaxOut4,xz3,xz4,xz5,colmut1,hx,k2,opfGoal,times,... secondBegin,Tn,Tn1,capKn,capKn1,kt0,q0,mut,mub,oneTrans,oneCap,para1,para2,para3,mismatch,ifconverge]=seting(nodeNum,lineNum,transNum,transK0,... capNum,capK,pvQmin,pvQmax,pvNum,pvNode,pvV,pgNum,pgNode,Pmin,Pmax,Pg,Qg,centrPara); disp (['迭代次数 互补间隙 最大不平衡量 最优潮流目标']) %% 进入循环体求解 for k=0:kmax [mut,mub,dNCPT,dNCPT_dt,d2NCPT_dt2,dT_Zt,transGmut,transBmut,dNCPcap,dNCP_dcap,d2NCP_dcap2,capBi_dBci,H1,N1,diagE,diagF,hx,Gap,max_hx,fx,... Pg,Qg,xz1,dmaxOut1,xz2,dmaxOut2,xz3,dmaxOut3,xz4,dmaxOut4,PGR,Tn,kt0,capKn,q0,Tn1,capKn1,T,G]=midResult(nodeNum,pgNode,opfGoal,z,l,u,w,Pg,... Qg,e,f,G,B,a,b,c,PG,Pd,Qd,k,xz1,dmaxOut1,xz2,dmaxOut2,xz3,dmaxOut3,xz4,dmaxOut4,times,mut,mub,oneTrans,oneCap,para3,accuracy,T,capK,transKmin,... transKmax,capKmin,capBi,mu,transNum,transG,transB,Tn,Tn1,capKn,capKn1,kt0,q0,mismatch); if Gapaccuracy; %判断靠近下界的离散量是否达到精度 dicVar1=[(Tn1-T);(capKn1-capK)]>accuracy; %判断靠近上界的离散量是否达到精度 ymut=dicVar.*dicVar1; if (sum(ymut)==0) %判断离散变量是否达到精度 ifconverge=1; break end end end mu=(centrPara*Gap)/(2*r); %计算扰动因子 if times==1 %预计算过程 %形成修正方程系数矩阵 [HSB,dhx,dgx,dgxT,N1,H1,diagE,diagF,GIJ,BIJ,L1,U1,LZ,UW]=formHSB(nodeNum,balNode,balNum,opfGoal,lineNum,transNum,lineI,lineJ,transI,... transJ,T,transG,transB,pgNode,pvNode,G,B,e,f,pgNum,pvNum,pgvNum,xNum,capNum,capI,a,m,r,l,u,z,w,y,H1,N1,diagE,diagF); %形成扰动KKT方程不平衡量矩阵 [const,Lz,Lw,L1,U1]=Constant(nodeNum,opfGoal,lineNum,transNum,transI,transJ,transG,transB,T,transKmin,transKmax,lineI,lineJ,balNode,... mu,G,B,GIJ,BIJ,e,f,Pmin,Pmax,a,b,PG,QR,pvQmin,pvQmax,pgNum,pgvNum,xNum,capNum,capK,capKmin,capKmax,pgNode,pvNode,dhx,dgx,L1,U1,LZ,... UW,l,u,y,r,capBi,Vmin,Vmax,hx); %求解方程并更新变量 [z,l,u,w,x,y,PG,QR,Pg,Qg,e,f,T,capK,G,B]=solve(HSB,const,mu,r,nodeNum,balNode,pgNum,pvNum,xNum,transNum,capNum,capI,capBi,capK,pgNode,... pvNode,m,x,y,z,l,u,w,Lz,Lw,LZ,UW,dgxT,L1,U1,G,B,transI,transJ,T,transG,transB,transG1,transB1,PG,QR,Pg,Qg); else %光滑处理过程 k2=k2+1; xz5=[xz5,k2]; %保留每次迭代光滑参数的值，用于作图 colmut1=[colmut1,mut(1)]; %形成修正方程系数矩阵 [HSB,dhx,dgx,dgxT,N1,H1,diagE,diagF,GIJ,BIJ,L1,U1,LZ,UW,transG,transB,mut,mub]=snd_formHSB(nodeNum,balNode,balNum,opfGoal,lineNum,... transNum,lineI,lineJ,transI,transJ,T,transG,transB,pgNode,pvNode,G,B,e,f,pgNum,pvNum,pgvNum,xNum,a,m,r,l,u,z,w,y,capNum,capI,... mut,mub,dNCPT_dt,d2NCPT_dt2,dT_Zt,transGmut,transBmut,dNCP_dcap,d2NCP_dcap2,H1,N1,diagE,diagF); %形成扰动KKT方程不平衡量矩阵 [const,Lz,Lw,L1,U1,hx,origConst,gxx,gsx]=snd_Constant(nodeNum,opfGoal,lineNum,transNum,transI,transJ,transG,transB,lineI,lineJ,... balNode,mu,G,B,GIJ,BIJ,e,f,Pmin,Pmax,a,b,PG,QR,pvQmin,pvQmax,pgNum,pgvNum,xNum,pgNode,pvNode,dhx,dgx,L1,U1,LZ,UW,l,u,y,r,... dT_Zt,capNum,z,w,Vmin,Vmax,hx); %求解方程并更新变量 [z,l,u,w,x,y,PG,QR,Pg,Qg,e,f,T,G,B,capK,meritF]=snd_solve(HSB,const,mu,r,nodeNum,balNode,pgNum,pvNum,pgvNum,xNum,transNum,pgNode,... pvNode,m,Lz,Lw,LZ,UW,dgxT,L1,U1,transI,transJ,transG,transB,transG1,transB1,Tn,capKn,capNum,capI,Tn1,capKn1,x,y,z,l,u,w,... G,B,PG,QR,Pg,Qg,opfGoal,Gama,Arlfa,meritF,k2,T,capK,lineNum,lineI,lineJ,Pd,Qd,gxx,gsx,mut,mub); end end if k>=kmax disp '最优潮流计算不收敛！' disp(k) else disp '最优潮流计算收敛！' disp(['迭代',num2str(k),'次！']) end toc; %% 输出计算结果 OutResult(G,B,balNode,e,f,T,k,Gap,fx,transI,transJ,filename,nodeNum,PGR,pgNode,QR,pvNode,capI,capK,firstT,firstcapK,Tn,Tn1,capKn,capKn1,... firstfx,toc,max_hx,ifconverge,Vfirst,Pgfirst,Qgfirst); %% 作互补间隙和最大不平衡量收敛曲线图 % xz=0:1:k; % plot(xz1,dmaxOut1,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3) % hold on; % plot(xz2,dmaxOut2,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3) % hold on; % xlabel('迭代次数') % ylabel('互补间隙') % plot(xz3,dmaxOut3,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3) % hold on; % plot(xz4,dmaxOut4,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3) % xlabel('迭代次数') % ylabel('最大不平衡量') %% 作mut在迭代过程中的变化图 % plot(xz5,colmut1,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3) % xlabel('第二阶段优化迭代次数') % ylabel('μ') % hold on; %% 作二次优化互补间隙收敛曲线 % plot(xz5,dmaxOut2,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3) % hold on; % xlabel('迭代次数') % ylabel('互补间隙')