基于内点法的最优潮流计算
- 格式:docx
- 大小:684.98 KB
- 文档页数:34
基于功率—电流混合潮流约束的内点法最优潮流江全元,黄志光(浙江大学电气工程学院,浙江省杭州市310027)摘要:提出了一种基于直角坐标下功率—电流混合型潮流约束的最优潮流模型。
该模型对系统中的非零注入功率节点采用功率失配型潮流约束,而对零注入功率节点采用电流失配型潮流约束。
这种混合模型结合了功率和电流型潮流方程的优点:对于零注入功率节点,该模型具有电流型潮流方程一阶导数为常数、二阶导数为0的特点,从而使雅可比矩阵和海森矩阵更容易计算;对于非零注入功率节点,该模型也比电流型潮流方程更好处理。
该模型特别适合非线性预测—校正内点法的最优潮流,多个大规模算例测试证明该模型收敛性更好,计算效率更高,尤其适合求解含较多零注入功率节点的大规模电力系统最优潮流问题。
关键词:最优潮流;内点法;电流失配;功率失配;混合模型中图分类号:TM744收稿日期:2008211219;修回日期:2009203202。
教育部科学技术研究重点项目(107063);浙江省自然科学基金资助项目(R1080089)。
0 引言自从20世纪60年代Carpentier 提出最优潮流(O PF )问题以来,该问题受到越来越多的关注[127]。
Dommel 和Tinney 建立了OPF 问题的非线性规划模型。
实际应用中,OPF 有不同的表述形式,学者们也提出了不同的算法求解该问题。
近年来,内点法及其改进算法由于其突出的计算性能在OPF 问题中得到了广泛的应用[8214]。
O PF 问题中,网络潮流约束可表示为直角坐标或极坐标形式。
文献[15]比较了直角坐标和极坐标2种形式的OPF 模型,指出二者有相似的收敛特性和计算效率。
但在O PF 的内点算法中,直角坐标系的表达更简洁[16]。
在潮流计算时,系统潮流方程通常表示为功率失配形式,而电流失配型潮流方程突出了电力网络本身是线性、而节点注入是非线性的特点[17],一般来说更加适合求解负荷潮流问题[18]。
基于Matlab符号计算工具箱的内点法最优潮流研究
李尹;韦化
【期刊名称】《电力自动化设备》
【年(卷),期】2003(023)007
【摘要】为提高最优潮流算法的通用性,利用Matlab符号计算工具箱完成了一种基于扰动KKT条件的内点算法最优潮流的符号计算.可以获得系统状态变量的显式符号结果,该方法使得复杂的最优潮流修正方程的形成与求解过程简化为在每次迭代中进行一次简单的代数替换.通过对4个不同的目标建模仿真,结果表明,该方法可极大地简化最优潮流计算程序的复杂程度,提高代码的通用性和易维护性.
【总页数】6页(P34-39)
【作者】李尹;韦化
【作者单位】广西大学,电气学院,广西,南宁,530004;广西大学,电气学院,广西,南宁,530004
【正文语种】中文
【中图分类】TM744;TP311
【相关文献】
1.基于内点法的最优潮流计算及算例分析 [J], 李春晓;何仁君
2.基于扰动KKT条件的原始-对偶内点法和分支定界法的最优潮流研究 [J], 范宏;韦化
3.基于加权预测-校正内点法的混合直流输电系统最优潮流 [J], 林子杰;黄为民;卫志农;孙国强;孙永辉
4.基于简化零空间内点法VSC-HVDC离散化最优潮流的研究 [J], 张昕;张勇;钱伟杰;章慧芸;胡仁
5.基于内点法的含风电场电力系统随机最优潮流 [J], 龚锦霞; 郑元黎
因版权原因,仅展示原文概要,查看原文内容请购买。
基于内点法的最优潮流计算GE GROUP system office room 【GEIHUA16H-GEIHUA GEIHUA8Q8-摘要内点法是一种能在可行域内部寻优的方法,即从初始内点出发,沿着中心路径方向在可行域内部直接走向最优解的方法。
其中路径跟踪法是目前最具有发展潜力的一类内点算法,该方法鲁棒性强,对初值的选择不敏感,在目前电力系统优化问题中得到了广泛的应用。
本文采用路径跟踪法进行最优求解,首先介绍了路径跟踪法的基本模型,并且结合具体算例,用编写的Matlab程序进行仿真分析,验证了该方法在最优潮流计算中的优越性能。
关键词:最优潮流、内点法、路径跟踪法、仿真目次0、引言........................................................1、路径跟踪法的基本数学模型....................................2、路径跟踪法的最优潮流求解思路................................3、具体算例及程序实现流程......................................3.1、算例描述..............................................3.2、程序具体实现流程......................................4、运行结果及分析..............................................4.1 运行结果..............................................4.2结果分析 ..............................................5、结论........................................................6、编程中遇到的问题............................................ 参考文献....................................................... 附录...........................................................0、引言电力系统最优潮流,简称OPF(Optimal Power Flow)。
基于ADMAT自动微分工具箱的内点法最优潮流计算鲍海波1,韦化1,2,莫东1(1.广西大学电力系统最优化研究所,2.广西电力系统最优化与节能技术重点实验室,广西壮族自治区南宁市530004)摘要:为了提高最优潮流程序的通用性和灵活性,避免计算内点法最优潮流时繁琐的导数公式推导,利用ADMAT自动微分工具箱完成了一种更为灵活的内点法最优潮流计算。
无需手动编程构造雅克比矩阵和海森矩阵,减少了编程出错可能性,且便于修改、增减最优潮流模型的目标函数和约束条件。
程序实现过程中,结合ADMAT工具箱的求解特点,研究了约束条件和目标函数的雅克比矩阵和海森矩阵的稀疏模式,对程序进一步优化。
多个系统的测试结果表明,自动微分技术应用于电力系统最优潮流计算的可行性和优越性,所设计的程序具有较高的计算效率。
关键词:自动微分;最优潮流;内点法;ADMAT0 引言经过四十多年的发展,最优潮流(optimal power flow ,OPF)技术在电力系统运行、规划、调度等领域得到了广泛的应用[1-2]。
现今,先后用于最优潮流计算方法有:牛顿法、内点法、遗传算法、差分进化算法等一系列方法[3-9],其中内点法具有多项式时间复杂性、计算效率高、鲁棒性好等优点,已经成功应用于最优潮流计算的工程实际。
内点法最优潮流计算过程中,需要手动推导目标函数和约束条件的1阶和2阶导数,以获得它们的雅克比(Jacobian)矩阵和海森(Hessian)矩阵。
这种手动编程方式具有很多弊端。
一方面,推导导数计算公式繁琐、容易出错,且代码不易于调试;另一方面,OPF在不同领域应用需要变化模型,增减或者修改部分约束、改变目标函数时,代码改动很大,限制了OPF程序代码的灵活性和可扩展性。
自动微分(automation differentiation,AD)技术的出现,成功解决了这个问题。
AD技术是机械的运用链式求导法则对计算机程序形式的函数求导的一种技术[10]。
基于内点法的最优潮流计算及算例分析
李春晓;何仁君
【期刊名称】《电气开关》
【年(卷),期】2018(055)001
【摘要】由于电力系统本身的复杂性,电力潮流优化具有规模大,约束条件多和非线性的特点.通过对最优潮流的求解,最终达到优化已有资源、降低发电厂耗量成本、减少电网线路损耗、提高电力系统输电能力等目标,其相比较传统的潮流计算具有良好的经济性.因此,最优潮流是电力系统中及受关注的课题,目前也有很多针对其做出的研究.本文综述了电力系统最优计算的数学模型和优化方法的研究现状.介绍了内点法的理论基石和基本原理,建立了最优潮流的数学模型,并对该模型采用内点法进行求解,最后通过实际算例加以验证.
【总页数】5页(P32-36)
【作者】李春晓;何仁君
【作者单位】广西大学电气工程学院,广西南宁 530004;广西大学电气工程学院,广西南宁 530004
【正文语种】中文
【中图分类】TM71
【相关文献】
1.基于Matlab符号计算工具箱的内点法最优潮流研究 [J], 李尹;韦化
2.基于改进内点法的含风电场的系统最优潮流计算 [J], 顾承红;艾芊
3.基于SolidWorks设计算例起吊系统的有限元分析 [J], 叶青玉
4.基于PSASP的矿区电网多目标最优潮流计算分析 [J], 唐翔; 黄倩
5.基于区间算法的配电网三相潮流计算及算例分析 [J], 王成山;王守相
因版权原因,仅展示原文概要,查看原文内容请购买。
内点法最优潮流MATLAB算法clear;%clc;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%数据加载n=input('请输入要计算的节点系统(5):')load Node5.txt;%节点数据load Branch5.txt;%支路数据load Generator5.txt;%发电机数据Node=Node5;Branch=Branch5;Generator=Generator5;%节点数据处理N=Node(:,1);%节点号Type=Node(:,2);%节点类型Uamp=Node(:,3);%节点电压幅值Dlta=Node(:,4);%节点电压相角Pd=Node(:,5);%节点负荷有功Qd=Node(:,6);%节点负荷无功Pg=Node(:,7);%节点出力有功Qg=Node(:,8);%节点出力无功Umax=Node(:,9);%节点电压幅值上限 Umin=Node(:,10);%节点电压幅值下限Bc=Node(:,11);%节点补偿电容电纳值 %支路数据处理Nbr=Branch(:,1);%支路号Nl=Branch(:,2);%支路首节点Nr=Branch(:,3);%支路末节点R=Branch(:,4);%支路电阻X=Branch(:,5);%支路电抗Z=R+1i*X;%支路阻抗=支路电阻+支路电抗 Bn=Branch(:,6);%支路对地电纳K=Branch(:,7);%支路变压器变比,0表示无变压器 Ptmax=Branch(:,8);%线路传输功率上限%发电机数据处理Ng=Generator(:,1);%发电机序号Nbus=Generator(:,2);%所在母线号Pumax=Generator(:,3);%发电机有功出力上界 Qumax=Generator(:,4);%发电机无功出力上界 Pumin=Generator(:,5);%发电机有功出力下界Qumin=Generator(:,6);%发电机无功出力下界a2=Generator(:,7);%燃料耗费曲线二次系数a1=Generator(:,8);%燃料耗费曲线一次系数a0=Generator(:,9);%燃料耗费曲线常数项%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%n=length(N);%节点个数ng=length(Ng);%发电机台数nbr=length(Nbr);%支路个数x=zeros(2*(ng+n),1);%控制变量+状态变量x(1:ng)=Pg(Nbus);x(ng+1:2*ng)=Qg(Nbus);x((2*ng+2):2:2*(ng+n))=Uamp; x((2*ng+1):2:2*(ng+n)-1)=Dlta; l=0.8*ones(2*ng+n+nbr,1);%松弛变量u=1.1*ones(2*ng+n+nbr,1);%松弛变量w=-1.5*ones(2*ng+n+nbr,1);%拉格朗日乘子z=ones(2*ng+n+nbr,1);%拉格朗日乘子y=zeros(2*n,1);%拉格朗日乘子y(1:2:2*n-1)=1e-3;y(2:2:2*n)=-1e-3;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算不等式约束的上下限%%%%%%%%%%%%%%%%%%%%%%%%%gmingmin=zeros(2*ng+n+nbr,1);gmin(1:ng)=Pumin;gmin(ng+1:2*ng)=Qumin;gmin(2*ng+1:2*ng+n)=Umin;gmin(2*ng+n+1:2*ng+n+nbr)=-Ptmax; %gmaxgmax=zeros(2*ng+n+nbr,1);gmax(1:ng)=Pumax;gmax(ng+1:2*ng)=Qumax;gmax(2*ng+1:2*ng+n)=Umax;gmax(2*ng+n+1:2*ng+n+nbr)=Ptmax;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成导纳矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Y=zeros(n,n);%%%%%%%%%%%%%%%%%%%%计算非对角元素%%%%%%%%%%%%%%%%%%%%% for ii=1:nbr if K(ii)==0%非变压器支路Y(Nl(ii),Nr(ii))=-1/Z(ii);Y(Nr(ii),Nl(ii))=Y(Nl(ii),Nr(ii));else%变压器支路Y(Nl(ii),Nr(ii))=-1/Z(ii)/K(ii);Y(Nr(ii),Nl(ii))= Y(Nl(ii),Nr(ii));endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算对角元素%%%%%%%%%%%%%%%%%%%%%%for ii=1:n%将支路导纳加入到对角元素中for jj=1:nbrif K(jj)==0&&(Nl(jj)==ii||Nr(jj)==ii)%非变压器支路Y(ii,ii)=Y(ii,ii)+1/Z(jj);else if K(jj)~=0&&(Nl(jj)==ii||Nr(jj)==ii)%变压器支路Y(ii,ii)=Y(ii,ii)+1/Z(jj)/K(jj);endendendendfor ii=1:nbr%将对地电纳加入到对角元素中if K(ii)==0%非变压器支路Y(Nl(ii),Nl(ii))=Y(Nl(ii),Nl(ii))+1i*Bn(ii);Y(Nr(ii),Nr(ii))=Y(Nr(ii),Nr(ii))+1i*Bn(ii);else%变压器支路Y(Nr(ii),Nr(ii))=Y(Nr(ii),Nr(ii))+(K(ii)-1)/K(ii)/Z(ii);Y(Nl(ii),Nl(ii))=Y(Nl(ii),Nl(ii))+(1-K(ii))/K(ii)/K(ii)/Z(ii);endendfor ii=1:nY(ii,ii)=Y(ii,ii)+i*Bc(ii);endG=real(Y);%电导B=imag(Y);%电纳%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%k=0;%迭代次数Kmax=150;%最大迭代次数iteration=1e-4;%误差精度delta=0.08;Gap=(l'*z-u'*w)*ones(Kmax,1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% while k<50%计算互补间隙GapGap(k+1)=l'*z-u'*w;if Gap>iterationmiu=delta*Gap(k+1)/(2*(2*ng+n+nbr)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%形成系数矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%相角差计算%%%%%%%%%%%%%%%%%%%%%%theta=zeros(n,n);for ii=1:nfor jj=1:ntheta(ii,jj)=Dlta(ii)-Dlta(jj);endend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%1、等式约束雅克比矩阵%%%%%%%%%%%%%%%%pxh=zeros(2*(ng+n),2*n); %%%%%%%%%%%%%%%%%%%%%%%ah/aP%%%%%%%%%%%%%%%%%%%%%%%for ii=1:ngpxh(Ng(ii),2*Nbus(ii)-1)=1;end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%ah/aQ%%%%%%%%%%%%%%%%%%%%%%%for ii=1:ngpxh(Ng(ii)+ng,2*Nbus(ii))=1;end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ah/ax%%%%%%%%%%%%%%%%%%%%%%%HH=zeros(n,n);JJ=zeros(n,n);NN=zeros(n,n);LL=zeros(n,n);for ii=1:nfor jj=1:nif ii~=jj%i!=j时的情况%非对角元素HH(ii,jj)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));JJ(ii,jj)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin (theta(ii,jj)));NN(ii,jj)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));LL(ii,jj)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角元素HH(ii,ii)=HH(ii,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));JJ(ii,ii)=JJ(ii,ii)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)) );NN(ii,ii)=NN(ii,ii)-Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));LL(ii,ii)=LL(ii,ii)-Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));endendNN(ii,ii)=NN(ii,ii)-2*Uamp(ii)*G(ii,ii);LL(ii,ii)=LL(ii,ii)+2*Uamp(ii)*B(ii,ii);endpxh(1+2*ng:2:2*(n+ng)-1,1:2:2*n-1)=HH';pxh(1+2*ng:2:2*(n+ng)-1,2:2:2*n)=JJ';pxh(2+2*ng:2:2*(n+ng),1:2:2*n-1)=NN';pxh(2+2*ng:2:2*(n+ng),2:2:2*n)=LL';%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2、不等式约束的雅克比矩阵%%%%%%%%%%%%%%%%%%%% %g1:电源有功出力上下限约束ag1aP=eye(ng,ng);ag1aQ=zeros(ng,ng);ag1ax=zeros(2*n,ng);%g2:电源无功出力上下限约束ag2aP=zeros(ng,ng);ag2aQ=eye(ng,ng);ag2ax=zeros(2*n,ng);%g3:节点电压幅值上下限约束ag3aP=zeros(ng,n);ag3aQ=zeros(ng,n);ag3ax=zeros(2*n,n);for ii=1:nag3ax(2*ii,ii)=1;end%g4:线路潮流上下限约束ag4aP=zeros(ng,nbr);ag4aQ=zeros(ng,nbr);ag4ax=zeros(2*n,nbr);for ii=1:nfor jj=1:nbrif Nl(jj)==iiag4ax(2*ii-1,jj)=-Uamp(Nl(jj))*Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*sin(theta(Nl(jj),N r(jj)))-B(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj))));ag4ax(2*ii,jj)=Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj )))+B(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj))))-2*Uamp(Nl(jj))*G(Nl(jj),Nr(jj));endif Nr(jj)==iiag4ax(2*ii-1,jj)=Uamp(Nl(jj))*Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr (jj)))-B(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj))));ag4ax(2*ii,jj)=Uamp(Nl(jj))*(G(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj )))+B(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj))));endendendpxg=[ag1aP ag2aP ag3aP ag4aP;ag1aQ ag2aQ ag3aQ ag4aQ;ag1ax ag2ax ag3ax ag4ax];%此即为不等式约束的雅克比矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%3、对角矩阵%%%%%%%%%%%%%%%%%%%%%%%% L_1Z=zeros(2*ng+n+nbr,2*ng+n+nbr);U_1W=zeros(2*ng+n+nbr,2*ng+n+nbr);for ii=1:2*ng+n+nbrL_1Z(ii,ii)=z(ii)/l(ii);U_1W(ii,ii)=w(ii)/u(ii);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%海森伯矩阵%%%%%%%%%%%%%%%%%%%%%%%%%% %将海森伯矩阵分为4块:H1,H2,H3,H4%%%%%%%%%%%%%%%%%%%%%H1%%%%%%%%%%%%%%%%%%%%%%A2=diag(a2);H1=zeros(2*(ng+n),2*(ng+n));H1(1:ng,1:ng)=2*A2;%%%%%%%%%%%%%%%%%%%%H2%%%%%%%%%%%%%%%%%%%%%%H2=zeros(2*(ng+n),2*(ng+n));A=zeros(2*n,2*n);Apb=zeros(2*n,2*n,n);Aqb=zeros(2*n,2*n,n);for ii=1:nfor jj=1:n %元素位置为:1 2if ii~=jj % 3 4%对角线上与ii对应的元素%ApApb(2*ii-1,2*ii-1,ii)=Apb(2*ii-1,2*ii-1,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(i i,jj)));%对角线处1号元素Apb(2*ii-1,2*ii,ii)=Apb(2*ii-1,2*ii,ii)+Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角线处2号元素%%3号元素与之相等%AqAqb(2*ii-1,2*ii-1,ii)=Aqb(2*ii-1,2*ii-1,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角线处1号元素Aqb(2*ii-1,2*ii,ii)=Aqb(2*ii-1,2*ii,ii)-Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%对角线处2号元素%%3号元素与之相等%对角线上与jj对应的元素%ApApb(2*jj-1,2*jj-1,ii)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(i i,jj)));%对角线处1号元素Apb(2*jj-1,2*jj,ii)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); %对角线处2号元素Apb(2*jj,2*jj-1,ii)=Apb(2*jj-1,2*jj,ii);%3号元素与2号元素相等%AqAqb(2*jj-1,2*jj-1,ii)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角线处1号元素Aqb(2*jj-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj ))); %对角线处2号元素Aqb(2*jj,2*jj-1,ii)=Aqb(2*jj-1,2*jj,ii);%3号元素与2号元素相等%4号元素为0%非对角线行元素%ApApb(2*ii-1,2*jj-1,ii)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)) );%非对角线行处1号元素Apb(2*ii-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处2号元素Apb(2*ii,2*jj-1,ii)=-Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处3号元素Apb(2*ii,2*jj,ii)=-(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处4号元素%AqAqb(2*ii-1,2*jj-1,ii)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处1号元素Aqb(2*ii-1,2*jj,ii)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处2号元素Aqb(2*ii,2*jj-1,ii)=Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处3号元素Aqb(2*ii,2*jj,ii)=-(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处4号元素%非对角线列元素%ApApb(2*jj-1,2*ii-1,ii)=Apb(2*ii-1,2*jj-1,ii);%非对角线列处1号元素Apb(2*jj-1,2*ii,ii)=Apb(2*ii,2*jj-1,ii);%非对角线列处2号元素Apb(2*jj,2*ii-1,ii)=Apb(2*ii-1,2*jj,ii);%非对角线列处3号元素Apb(2*jj,2*ii,ii)=Apb(2*ii,2*jj,ii);%%非对角线列处4号元素%AqAqb(2*jj-1,2*ii-1,ii)=Aqb(2*ii-1,2*jj-1,ii);%非对角线列处1号元素Aqb(2*jj-1,2*ii,ii)=Aqb(2*ii,2*jj-1,ii);%非对角线列处2号元素Aqb(2*jj,2*ii-1,ii)=Aqb(2*ii-1,2*jj,ii);%非对角线列处3号元素Aqb(2*jj,2*ii,ii)=Aqb(2*ii,2*jj,ii);%%非对角线列处4号元素endend%对角线上与ii对应的元素%ApApb(2*ii,2*ii-1,ii)=Apb(2*ii-1,2*ii,ii);%对角线处3号元素与2号元素相等Apb(2*ii,2*ii,ii)=-2*G(ii,ii);%对角线处4号元素%AqAqb(2*ii,2*ii-1,ii)=Aqb(2*ii-1,2*ii,ii);%对角线处3号元素与2号元素相等Aqb(2*ii,2*ii,ii)=2*B(ii,ii);%对角线处4号元素endfor ii=1:nA=A+Apb(:,:,ii)*y(2*ii-1)+Aqb(:,:,ii)*y(2*ii);endH2(2*ng+1:2*(ng+n),2*ng+1:2*(ng+n))=A;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%H3%%%%%%%%%%%%%%%%%%%%%%H3=zeros(2*(ng+n),2*(ng+n));A3=zeros(2*n,2*n);Apc=zeros(2*n,2*n,nbr);for ii=1:nbr%对角线上iiApc(2*Nl(ii)-1,2*Nl(ii)-1,ii)=-Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))+B( Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii))));Apc(2*Nl(ii)-1,2*Nl(ii),ii)=-Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii))));Apc(2*Nl(ii),2*Nl(ii)-1,ii)=Apc(2*Nl(ii)-1,2*Nl(ii),ii);Apc(2*Nl(ii),2*Nl(ii),ii)=-2*G(Nl(ii),Nr(ii));%对角线上jjApc(2*Nr(ii)-1,2*Nr(ii)-1,ii)=-Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))+B( Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii))));Apc(2*Nr(ii)-1,2*Nr(ii),ii)=Uamp(Nl(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii))));Apc(2*Nr(ii),2*Nr(ii)-1,ii)=Apc(2*Nr(ii)-1,2*Nr(ii),ii);Apc(2*Nr(ii),2*Nr(ii),ii)=0;%非对角线ijApc(2*Nl(ii)-1,2*Nr(ii)-1,ii)=Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii )))+B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii))));Apc(2*Nl(ii)-1,2*Nr(ii),ii)=-Uamp(Nl(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii))));Apc(2*Nl(ii),2*Nr(ii)-1,ii)=Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii))));Apc(2*Nl(ii),2*Nr(ii),ii)=G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii))) +B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)));%非对角线jiApc(2*Nr(ii)-1,2*Nl(ii)-1,ii)=Apc(2*Nl(ii)-1,2*Nr(ii)-1,ii);Apc(2*Nr(ii)-1,2*Nl(ii),ii)=Apc(2*Nl(ii),2*Nr(ii)-1,ii);Apc(2*Nr(ii),2*Nl(ii)-1,ii)=Apc(2*Nl(ii)-1,2*Nr(ii),ii);Apc(2*Nr(ii),2*Nl(ii),ii)=Apc(2*Nl(ii),2*Nr(ii),ii);%求和c=z+w;A3=A3+Apc(:,:,ii)*c(2*ng+n+ii);endH3(2*ng+1:2*(ng+n),2*ng+1:2*(ng+n))=A3;%%%%%%%%%%%%%%%%%%%%%%%H4%%%%%%%%%%%%%%%%%%%%%%%%%H4=pxg*(L_1Z-U_1W)*pxg';%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%H=-H1+H2+H3-H4;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成常数项%%%%%%%%%%%%%%%%%%%%%%%%% %Lyh=zeros(2*n,1);for ii=1:nh(2*ii-1)=Pg(ii)-Pd(ii);h(2*ii)=Qg(ii)-Qd(ii);for jj=1:nh(2*ii-1)=h(2*ii-1)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(t heta(ii,jj)));h(2*ii)=h(2*ii)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));endendLy=h;%Lz%g(x)gx=zeros(2*ng+n+nbr,1);gx(1:ng)=x(1:ng);gx(ng+1:2*ng)=x(ng+1:2*ng);gx(2*ng+1:2*ng+n)=x(2*ng+2:2:2*(ng+n));for ii=1:nbrgx(2*ng+n+ii)=Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta( Nl(ii),Nr(ii)))+B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii))))-Uamp(Nl(ii))*Uamp(Nl(ii))*G(Nl(ii),Nr(ii));endLz=gx-l-gmin;%LwLw=gx+u-gmax;%Lle=ones(2*ng+n+nbr,1);LZ=zeros(2*ng+n+nbr,2*ng+n+nbr);for ii=1:2*ng+n+nbr;LZ(ii,ii)=l(ii)*z(ii);endLl=LZ*e-miu*e;%LuUW=zeros(2*ng+n+nbr,2*ng+n+nbr);for ii=1:2*ng+n+nbrUW(ii,ii)=u(ii)*w(ii);endLu=UW*e+miu*e;%Lx'Lx1=zeros(2*(ng+n),1);Lx1(1:ng)=2*a2.*x(1:ng)+a1;Lx2=pxh*y;Lx3=pxg*c;Lx41=zeros(2*(ng+n),1);Lx42=zeros(2*(ng+n),1);for ii=1:2*ng+n+nbrLx41(ii)=(Ll(ii)+z(ii)*Lz(ii))/l(ii);Lx42(ii)=(Lu(ii)-w(ii)*Lw(ii))/u(ii);endLx4=pxg*(Lx41+Lx42);Lx=Lx1-Lx2-Lx3;Lxx=Lx+Lx4; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%求出修正量%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %dx,dyHxy=[H pxh;pxh' zeros(2*n,2*n)];LxLy=[Lxx;-Ly];I=eye(2*(ng+n)+2*n);dxdy=I/Hxy*LxLy;dx=dxdy(1:2*(ng+n));dy=dxdy(2*(ng+n)+1:2*(ng+n)+2*n);%dldl=pxg'*dx+Lz;%dudu=-pxg'*dx-Lw;%dzdz=zeros(2*ng+n+nbr,1);for ii=1:2*ng+n+nbrdz(ii)=(-Ll(ii)-z(ii)*dl(ii))/l(ii);end%dwdw=zeros(2*ng+n+nbr,1);for ii=1:2*ng+n+nbrdw(ii)=(-Lu(ii)-w(ii)*du(ii))/u(ii);end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%计算alfap和alfad%%%%%%%%%%%%%%%%%%%%%%%% alfap=1;alfad=1;for ii=1:2*ng+n+nbrif dl(ii)<0&&-l(ii)/dl(ii)<alfapalfap=-l(ii)/dl(ii);endif du(ii)<0&&-u(ii)/du(ii)<alfapalfap=-u(ii)/du(ii);endif dz(ii)<0&&-z(ii)/dz(ii)<alfadalfad=-z(ii)/dz(ii);endif dw(ii)>0&&-w(ii)/dw(ii)<alfadalfad=-w(ii)/dw(ii);endendalfap=0.9995*alfap;alfad=0.9995*alfad;x=x+alfap*dx;l=l+alfap*dl;u=u+alfap*du;y=y+alfad*dy;z=z+alfad*dz;w=w+alfad*dw;%迭代功率、电压幅值和相角for ii=1:ngPg(Nbus(ii))=x(ii);Qg(Nbus(ii))=x(ng+ii);endfor ii=1:nUamp(ii)=x(2*(ng+ii));Dlta(ii)=x(2*(ng+ii)-1);end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% k=k+1;elsebreak;endendfcost=0;for ii=1:ngfcost=fcost+a2(ii)*Pg(Nbus(ii))*Pg(Nbus(ii))+a1(ii)*Pg(Nbus(ii))+a0( ii);endfcostkplot(0:k,Gap(1:k+1),':*');PgQgUampfor ii=1:nif Type(ii)==3Dlta=Dlta-Dlta(ii)*ones(n,1);endendDlta。
基于内点法的最优潮流计算Company Document number:WTUT-WT88Y-W8BBGB-BWYTT-19998摘要内点法是一种能在可行域内部寻优的方法,即从初始内点出发,沿着中心路径方向在可行域内部直接走向最优解的方法。
其中路径跟踪法是目前最具有发展潜力的一类内点算法,该方法鲁棒性强,对初值的选择不敏感,在目前电力系统优化问题中得到了广泛的应用。
本文采用路径跟踪法进行最优求解,首先介绍了路径跟踪法的基本模型,并且结合具体算例,用编写的Matlab程序进行仿真分析,验证了该方法在最优潮流计算中的优越性能。
关键词:最优潮流、内点法、路径跟踪法、仿真目次0、引言电力系统最优潮流,简称OPF(Optimal Power Flow)。
OPF问题是一个复杂的非线性规划问题,要求满足待定的电力系统运行和安全约束条件下,通过调整系统中可利用控制手段实现预定目标最优的系统稳定运行状态。
针对不同的应用,OPF模型课以选择不同的控制变量、状态变量集合,不同的目标函数,以及不同的约束条件,其数学模型可描述为确定一组最优控制变量u,以使目标函数取极小值,并且满足如下等式和不等式。
{min u f(x,u)S.t.ℎ(x,u)=0g(x,u)≤0(0-1)其中min u f(x,u)为优化的目标函数,可以表示系统运行成本最小、或者系统运行网损最小。
S.t.ℎ(x,u)=0为等式约束,表示满足系统稳定运行的功率平衡。
g(x,u)≤0为不等式约束,表示电源有功出力的上下界约束、节点电压上下线约束、线路传输功率上下线约束等等。
电力系统最优潮流算法大致可以分为两类:经典算法和智能算法。
其中经典算法主要是指以简化梯度法、牛顿法、内点法和解耦法为代表的基于线性规划和非线性规划以及解耦原则的算法,是研究最多的最优潮流算法, 这类算法的特点是以一阶或二阶梯度作为寻找最优解的主要信息。
智能算法主要是指遗传算法和模拟退火发等,这类算法的特点是不以梯度作为寻优信息,属于非导数的优化方法。
因此经典算法的优点是能按目标函数的导数信息确定搜索方向,计算速度快,算法比较成熟,结果可信度高。
缺点是对目标函数及约束条件有一定的限制,可能出现局部极小时难以收敛。
而智能算法的优点是计算与导数无关,灵活性高,随机性强,缺点是算法不稳定,结果不可信,并且控制参数需凭经验给出。
通过对这些常见算法的简单比较,内点法具有其优越的性能,特别是路径跟踪法,其算法收敛迅速,鲁棒性强,对初值的选择不敏感,其迭代次数与系统规模或控制变量的数目关系不大,因此本文采用该方法进行最优计算。
1、路径跟踪法的基本数学模型内点法最初的基本思路是希望通过寻优迭代过程始终在可行域内进行,因此,初始点应在可行域内,并在可行域的边界设置‘障碍’使迭代点接近边界时其目标函数迅速增大,从而保证迭代点均在可行域的内点。
但是对于大规模实际问题而言,寻找初始点往往十分困难。
为此许多学者长期以来致力于内点算法初始‘内点’条件的改进。
以下介绍的路径跟踪法只要求在寻优过程中松弛变量和拉格朗日乘子满足简单的大于0或者小于0的条件,即可代替原来必须在可行域内求解的要求,使得计算过程大为简化。
一般可以将最优潮流模型简化为如下的非线性优化模型。
Obj. min u f(x,u)(1-1). S.t.ℎ(x,u)=0(1-2)g−≤g(x,u)≤g−(1-3)其中min u f(x,u)为优化的目标函数, S.t.ℎ(x,u)=0为等式约束, g(x,u)为不等式约束,路径跟踪内点法的基本思路是:首先将式(1-3)的不等约束变成等式约束:g (x,u )+u =g − (1-4) g (x,u )−l =g − (1-5)其中松弛变量 l =[l 1,…,l r ]T, u =[u 1,…,u r ]T,应满足u>0,l>0这样原问题就转化为问题A :Obj. min u f (x,u ). ()0()()h x g x u g g x l g=+=-=然后,把目标函数改造成障碍函数,该函数在可行域内应接近于原函数f(x),而在边界时变得很大。
一次可得带优化问题B :obj. 11min.()log()log()rrr r j j f x u l u u ==--∑∑. ()0()()h x g x u g g x l g=+=-=其中扰动因子或者障碍因子u>0。
当l 或u 接近边界时,以上函数将趋于无穷大,因此满足以上障碍目标函数的极小解不可能在边界上找到。
这样就通过目标函数的变化把含不等式限制的优化问题A 变成只含等式限制优化的问题B 了,因此可以直接用拉格朗日乘子法来求解。
优化问题B 的拉格朗日目标函数为:11()()[()][()]log()log()r rTTTr r j j L f x y h x z g x l g w g x u g u l u u ===-----+---∑∑ 式中:y ,z 和w 均为拉格朗日乘子。
因此最后简化的求解问题就是求取上述表达式的极小解。
2、 路径跟踪法的最优潮流求解思路路径跟踪法的最优潮流求解过程就是对拉格朗日目标函数求极小值问题: 式中:y ,z 和w 均为拉格朗日乘子。
该问题极小值存在的必要条件是拉格朗日函数对所有变量及乘子的偏导数为0。
即:11()()()(+w)0()0()0()+u 00U 0μμμμμμ--∂==∇-∇-∇=∂∂===∂∂==--=∂∂==-=∂∂==-⇒=-=∂∂==--⇒=+=∂x x x x y z w l lu uLL f x h x y g x z x L L h x y L L g x l g z L L g x g w L L z L e L LZe e l L L w e L UWe e u (2-1)通过上述表达式可以解出:2T T l z u w rμ-= (2-2)定义:T T Gap l z u w =-,称为互补间隙。
可得:2Gaprμ=(2-3)如果x *是优化问题A 的最优解,当u 固定时,x(u)是优化问题B 的解,那么当Gap →0,u →0时,产生的序列{x(u)}收敛至x *。
建议采用:2Gaprμσ=。
式中(0,1)σ∈称为中心参数,一般取,在大多数场合可获得较好的收敛效果。
通过偏导数为0的表达式可以可得内点法的修正方程为:111'0000000()000000000()00000'()0000()0l T z x u Tw x x x T y x L L z I L ZL l I g x U L w I L u I g x L x H h x L y h x μμ----∆⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥∆-∇⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-∆⎢⎥=⎢⎥⎢⎥⎢⎥-∆∇⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥∆∇⎢⎥⎢⎥⎢⎥-∆∇⎢⎥⎣⎦⎣⎦⎣⎦(2-4) 求解方程可得到第k 次迭代的修正量,于是最优解的一个新的近似解为:(1)()(1)()(1)()(1)()(1)()(1)()k k p k k p k k p k k d k k d k k d x x a x l l a l u u a u yya yz z a z w w a w++++++=+∆=+∆=+∆=+∆=+∆=+∆ (2-5)式中,p a 和d a 为步长:0.9995min{min(,0;,0),1}i i p i i i il ua l u l u --=∆<∆<∆∆ 0.9995min{min(,0;,0),1}i i d i i i iz wa z w z w --=∆<∆>∆∆ (2-6) 其潮流计算的流程图如下图1所示,其中初始化部分包括: (1)、设置松弛变量l 和u ,保证[l ,u]T >0(2)、设置拉格朗日乘子w 、y 、z ,满足[w<0,z>0,Y !=0]T (3)、设置优化问题的初值。
(4)、取中心参数(0,1)σ∈,给定计算精度,迭代次数初值K=0。
3、 具体算例及程序实现流程这部分主要有算例描述以及程序的实现流程两部分,其中算例描述主要是对系统参数以及优化问题进行说明。
而程序的实现流程主要描述的是最优潮流计算中所涉及的矩阵方程的描述。
、算例描述该算例为王锡凡编写的《现代电力系统分析》中的3-1的例题,是以系统燃料最省为最优潮流的目标函数。
选择该题目作为算例分析的原因是,该题目有比较详细的解题思路以及列写出了比较详细的迭代结果,方便对编写程序的运行结果进行比对。
求如下图所示简化系统的系统燃料最省的最优潮流计算:除了由上图所提供的系统母线负荷功率数据、线路参数和变压器之路参数数据、变压器便比数据之外,以下顺序给出了线路传输功率边界(表3-1),发电机有功无功出力上下界和燃料耗费曲线 参数(表3-2)。
若不作特殊说明,2435所有数据都是以标幺值形式给出,功率基准值为100MVA,母线电压上下界分别为和。
表3-1线路传输功率边界表3-2 发电机数据、程序具体实现过程针对上述系统,首先我们先列写出该算例的数学模型和有关计算公式。
在该算例中,共有5个节点,相应的状态量为:系统中有2台发电机,没有其他无功源,因此控制变量为:应该指出,此处发电机和无功源的编号与及诶单编号无关,是独立编号的。
这是因为系统中一个节点可能接有多台发电机的缘故。
因此系统中总变量共有14个:最优潮流的数学模型为:目标函数:约束条件:每个节点有两个潮流方程,共有10个等式约束条件,对非发电机而言:5151(cos sin )0(sin cos )0θθθθ===--+==-+-=∑∑i Di i j ij ij ij ij j i Di i j ij ij ij ij j P P V V G B Q Q V V G B (i=1,2,3)对发电机节点:5151(cos sin )0(sin cos )0θθθθ∈=∈==--+==-+-=∑∑∑∑i Gk Di i j ij ij ij ij k ij i Gk Di i j ij ij ij ij k ij P P P V V G B Q Q Q V V G B (i=4,5)式中:∈k i 表示第k 台发电机接在节点i 上。
不等式约束共有14个条件,分别是:根据以上模型可以形成修正方程。
该方程包括形成等式左边的系数矩阵和等式右边的常数项两部分。
1、形成系数矩阵1)、等式约束的雅克比矩阵 等式右端包括3个子矩阵: 其中: 其中:式中:i 为发电机的序号;j 为节点号;()∈i j 表示第i 台发电机是在节点j 上的。
111111111111111155551111555511111010.................θθθθθθθθ⨯∂∆∂∆∂∆∂∆⎡⎤⎢⎥∂∂∂∂⎢⎥⎢⎥∂∆∂∆∂∆∂∆⎢⎥∂∂∂∂⎢⎥∂⎢⎥=∂⎢⎥∂∆∂∆∂∆∂∆⎢⎥⎢⎥∂∂∂∂⎢⎥∂∆∂∆∂∆∂∆⎢⎥⎢⎥∂∂∂∂⎣⎦P Q P Q P P P P V V V V h x P Q P Q P Q P Q V V V V (潮流计算中的雅克比矩阵)2)、不等式约束的雅克比矩阵式中:1g 、2g 、3g 和4g 依次表示电源有功出力的上下界约束,无功电源出力的上下界约束,节点电压赋值的上下界约束和线路潮流约束。