基于内点法最优潮流计算 PPT
基于功率—电流混合潮流约束的内点法最优潮流

基于功率—电流混合潮流约束的内点法最优潮流江全元,黄志光(浙江大学电气工程学院,浙江省杭州市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]。
一种基于Karmarkar内点法的最优潮流算法

一种基于Karmarkar内点法的最优潮流算法
郝玉国;刘广一;于尔铿
【期刊名称】《中国电机工程学报》
【年(卷),期】1996(16)6
【摘要】以原-对偶内点算法(Karmarkar内点法的一种变形)为基本算法解算最优潮流问题,综合考虑非线性目标函数和约束条件,结合牛顿法最优潮流先进的稀疏矩阵技术,并且提出了一种新的原-对偶内点算法迭代步长选取原则和障碍参数修正策略。
算例表明本算法有较好的数值稳定性,优化结果精确,对不等式约束有较强的处理能力,显示了内点算法应用于大规模电力系统优化问题的良好前景。
【总页数】4页(P409-412)
【关键词】内点算法;最优潮流;电力系统
【作者】郝玉国;刘广一;于尔铿
【作者单位】电力科学研究院电网所
【正文语种】中文
【中图分类】TM744
【相关文献】
1.基于改进多中心校正解耦内点法的动态最优潮流并行算法 [J], 简金宝;杨林峰;全然
2.基于非线性多中心校正内点法的最优潮流算法 [J], 蔡广林;张勇军;任震
3.基于自动微分技术的内点法最优潮流算法 [J], 耿光超;江全元
4.基于退火粒子群和内点法的改进最优潮流算法 [J], 陈丽光;文波;聂一雄
5.基于预测校正内点法的最优潮流算法 [J], 陆子强
因版权原因,仅展示原文概要,查看原文内容请购买。
基于内点法最优潮流计算.32页PPT

36、如果我们国家的法律中只有某种 神灵, 而不是 殚精竭 虑将神 灵揉进 宪法, 总体上 来说, 法律就 会更好 。—— 马克·吐 温 37、纲纪废弃之日,便是暴政兴起之 时。— —威·皮 物特
38、若是没有公众舆论的支持,法律 是丝毫 没有力 量的。 ——菲 力普斯 39、一个判例造出另一个判例,它们 迅速累 聚,进 而变成 法律—爱献 生
6、最大的骄傲于最大的自卑都表示心灵的最软弱无力。——斯宾诺莎 7、自知之明是最难得的知识。——西班牙 8、勇气通往天堂,怯懦通往地狱。——塞内加 9、有时候读书是一种巧妙地避开思考的方法。——赫尔普斯 10、阅读一切好书如同和过去最杰出的人谈话。——笛卡儿
Thank you
潮流的计算机算法ppt课件

2020/3/31
1
带有最优 乘子的牛 顿潮流算
法
牛拉 法
保留非线 性直角坐
标法
保留非线 性直角坐 标快速潮
流法
简
简化
化
满足初 始条件 时为等 效算法
PQ分 解法
定雅克 比牛顿
法
2020/3/31
基本潮流
最优潮流牛顿算法
最优潮流简化梯度 算法
优化潮流
2
算法名称
算法特性
最优乘子法
能够有效地解决病态系统的潮流计算,且 永远不换发散
2020/3/31
3
内容提要
功率方程 牛拉法 P-Q分解法 保留非线性潮流算法 最小化潮流算法 最优潮流 潮流计算中稀疏技术的运用
2020/3/31
4
➢功率方程
电力系统中已知的往往是功率,需要用已知的功率来代替未
知的电流:
S%i
Pi
代入H、L的表达式
i j时
H ij
Pi
j
UiU j (Gij sin ij
Bij cosij )
cosij 1,Gij sinij 0 U iU j Bij
Lij
Pi U j
Uj
UiU j (Gij
sin ij
Bij
cosij )
cosij 1,Gij sinij 0 U iU j Bij
2020/3/31
15
i j时
H ii
Pi
i
Qi
U
2 i
Bii
U
B2
i ii
Lii
Qi
U
2 i
Bii
U
2 i
内点法最优潮流MATLAB算法

内点法最优潮流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。
潮流计算上机指导ppt课件

创建开发项目
打开“开始->程 序 ->Microsoft Visual Studio2008 -> Microsoft Visual Studio2008”,启动 Microsoft Visual Studio2008的主程序 界面如右
点击“文 件->新建>项目”, 启动了创 建项目对 话框如右
1.2.1 头文件和命名空间的引用 1.2.2 iostream类的使用
头文件和命名空间的引用
根据C语言标准,所有类和函数都是使用 头文件进行定义和说明的。在程序的开始 需要加入必需的头文件(.h)。C++类库 还增加了命名空间(namespace),程序所 用到的大多数类都包含在 “std”命名空间中,这些 都 需要在程序中加以说明。
成绩
提交课程设计论文报告,包含三部分内
容:程序设计体会、思考题答案和潮流程序 及计算结果。
比试。
1 基础知识
开发潮流计算程序使用的平Microsoft Visual Studio系列开发软件,本课程主要 采用Visual C/C++软件开发。下面介绍一 些常用的知识。 1.1创建开发项目 1.2简单输入输出程序测试 1.3程序的调试方法 1.4文件输入输出程序测试 1.5解方程程序测试
课程安排
课程共分为几个阶段: 程序设计初步知识(2学时) 数据结构、主程序框架、数据录入(2学时) 节点导纳矩阵的生成(2学时) 设置电压初始值、计算不平衡量(2学时) 建立雅可比矩阵、解修正方程(2学时) 设定收敛条件,完成潮流迭代程序(2学时) 计算线路潮流,网损,平衡节点功率,输出计算 结果(2学时)
// TEST.cpp: 主项目文件。 #include "stdafx.h" #include "fstream" #include "iostream" #include "NEquation.h" using namespace System; using namespace std; int main(array<System::String ^> ^args) { NEquation ObNEquation; //建立一个对象 ObNEquation.SetSize(2); //设置矩阵的维数为2维 ifstream infile1; infile1.open("Xishu.txt"); for(int i=0;i<2;i++) for(int j=0;j<2;j++) infile1>>ObNEquation.Data(i,j); infile1.close(); ifstream infile2; infile2.open("Hanshuzhi.txt"); for(int i=0;i<2;i++) infile2>>ObNEquation.Value(i); infile2.close(); ObNEquation.Run(); //调用NEquation头文件中的Run函数。 for(int i=0;i<2;i++) cout<<ObNEquation.Value(i)<<endl; //输出结果 Console::WriteLine(L"Hello World"); return 0; }
基于ADMAT 自动微分工具箱的内点法最优潮流计算—第27届高校电力系统及其自动化专业年会

基于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]。
基于内点法的最优潮流计算

基于内点法的最优潮流计算摘要内点法是一种能在可行域内部寻优的方法,即从初始内点出发,沿着中心路径方向在可行域内部直接走向最优解的方法。
其中路径跟踪法是目前最具有发展潜力的一类内点算法,该方法鲁棒性强,对初值的选择不敏感,在目前电力系统优化问题中得到了广泛的应用。
本文采用路径跟踪法进行最优求解,首先介绍了路径跟踪法的基本模型,并且结合具体算例,用编写的Matlab程序进行仿真分析,验证了该方法在最优潮流计算中的优越性能。
关键词:最优潮流、内点法、路径跟踪法、仿真目次0、引言 (1)1、路径跟踪法的基本数学模型 (2)2、路径跟踪法的最优潮流求解思路 (3)3、具体算例及程序实现流程 (7)3.1、算例描述 (7)3.2、程序具体实现流程 (8)4、运行结果及分析 (12)4.1 运行结果 (12)4.2结果分析 (18)5、结论 (19)6、编程中遇到的问题 (20)参考文献 (22)附录 (23)0、引言电力系统最优潮流,简称OPF(Optimal Power Flow)。
OPF问题是一个复杂的非线性规划问题,要求满足待定的电力系统运行和安全约束条件下,通过调整系统中可利用控制手段实现预定目标最优的系统稳定运行状态。
针对不同的应用,OPF模型课以选择不同的控制变量、状态变量集合,不同的目标函数,以及不同的约束条件,其数学模型可描述为确定一组最优控制变量u,以使目标函数取极小值,并且满足如下等式和不等式。
(0-1)其中为优化的目标函数,可以表示系统运行成本最小、或者系统运行网损最小。
为等式约束,表示满足系统稳定运行的功率平衡。
为不等式约束,表示电源有功出力的上下界约束、节点电压上下线约束、线路传输功率上下线约束等等。
电力系统最优潮流算法大致可以分为两类:经典算法和智能算法。
其中经典算法主要是指以简化梯度法、牛顿法、内点法和解耦法为代表的基于线性规划和非线性规划以及解耦原则的算法,是研究最多的最优潮流算法, 这类算法的特点是以一阶或二阶梯度作为寻找最优解的主要信息。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
8300
8200
8100
8000
7900
7800
7ห้องสมุดไป่ตู้00
7600
0
2
4
6
8 10 12 14 16
迭代次数
5节点目标函数变化曲线
102
0
10
10-2
10-4
-6
10
10-8
-10
10
0
2
4
6
8
10
12
14
16
迭代次数
5节点最大不平衡量变化曲线
目标函数
最大不平衡量
1092
1091.5
1091
1090.5
r
r
L f( x ) y T h ( x ) z T [ g ( x ) l g ] w T [ g ( x ) u g ] ulo lr ) u gl( o u r )g(
j 1
j 1
用牛顿法求解KKT方程,得到最优解。
L 0 , L 0 , L 0 , L 0 , L 0 , L 0 x y z w l u
1:1.05 2 0.08+0.30j 4 0.015j
1.05:1
3 0.03j
5
2+1j
j0.25
0.04+0.25j 0.25j
j0.25 3.7+1.3j
0.1+0.35j
1
1.6+0.8j
1+0.35j
2
7
0.0625j 8 0.0085+0.072j
0.0119+0.1008j 6 0.0586j
3
0.153j 0.032+0.161j
0.0745j 0.1045j
0.179j 0.039+0.017j
1.25+0.5j
9
0.088j 0.01+0.085j
5
0.9+0.3j
0.017+0.092j
0.079j 4 0.0576j
1
5节点系统结构图
9节点系统结构图
1、模型
5节点算例求解过程
内点法的优越性:
• 1、收敛速度快。 • 2、对系统规模不敏感。 • 3、对初始点不敏感。
二、最优潮流的原对偶内点算法
数学模型:
obj . min . f ( x ) s.t . h ( x ) 0
g g(x) g
f(x)为目标函数;h(x)为等式约束条件;g(x)为不等式约束条件。
原对偶内点算法:
定义对偶间隙和障碍参数为: GaplTzuTw
u Gap
2r
内点法小结
• 内点法实质上是牛顿法、对数壁垒函数法以及拉格朗日函 数法三者的结合。用对数壁垒函数处理不等式约束,用拉 格朗日函数处理等式约束,用牛顿法求解修正方程。
• (1)初始点的选取:跟踪中心轨迹内点法对初始点无要 求。
• (2)迭代收敛判据:对偶间隙小于某一给定值(最大潮 流偏差小于某一给定值)。
1090
1089.5
1089
1
2
3
4
5
6
7
8
9
10
迭代次数
9节点目标函数变化曲线
目标函数
4
3.9
3.8
3.7
3.6
3.5
3.4
3.3
3.2
3.1
3
2
4
6
8
10
12
迭代次数
30节点目标函数变化曲线
0
0
10
10
-2
10
10-4
10-6 1 2 3 4 5 6 7 8 9 10 11 迭代次数
9节点最大不平衡量变化曲线
最大不平衡量
-2
10
10-4
-6
10
0
2
4
6
8
10
12
迭代次数
30节点最大不平衡量变化曲线
最大不平衡量
收敛特性分析
2
2
2
10
10
10
100
100
100
-2
-2
-2
10
10
10
Gap Gap Gap
10-4
10-4
10-4
-6
-6
-6
10
10
10
10-8
0
2
4
6
8 10 12 14 16
迭代次数
5节点系统对偶间隙变化曲线
收敛特性分析
9节点系统迭代步长 30节点系统迭代步长
收敛特性分析
下表为计算过程中5节点系统的迭代步长:
分析可知,迭代过程开始时步长过小是制约5节点系统求解速 率的主要原因。
仿真结果分析
运用powerworld仿真的5节点算例结果如下图所示:
1.07 pu
1.10 pu
0.48 rad
0.40 rad
大家应该也有点累了,稍作休息
大家有疑问的,可以询问和交流
算法流程图:
初始化
计算互补间隙Gap 是
Gap<
否 计算扰动因子miu
求解修正方程,得各修正量△x,△y,△l,△v,△z,△w
计算步长ap和ad
更新原始变量和对偶变量
是 否
k<50
输出“计算不收敛!”
输出最 优解。
算例结构图
运用MATLAB最优潮流内点算法程序测试的5节点、9节点(30节点)d 等系统的结构图如下所示。
基于内点法最优潮流计算
主要内容
1、 课题研究的意义和现状 2、 最优潮流的原对偶内点算法 3、 最优潮流的预测校正内点算法 4、 结论
一、课题研究的意义和现状
概念:
最优潮流问题(OPF)就是在系统结构参数及负荷 给定的情况下,通过优选控制变量,确定能满足所 有的指定约束条件,并使系统的某个性能指标达到 最优时的潮流分布。
5节点算例求解过程
2、形成系数矩阵
5节点算例求解过程
5节点算例求解过程
3、形成常数项
算例迭代过程分析
PG 1
PG 2
PG 3
QG1
QG 2
QG3
各有功、无功电源出力随迭代次数的变化情况
算例迭代过程分析
1
V1
2
V2
3
V3
节点电压相角、幅值随迭代次数的变化情况
收敛特性分析
目标函数
8400
意义:
电力系统的经济运行一直是研究者们的热门课题。 随着人们对电能质量和安全性问题的重视,迫切需 要将三方面的要求统一起来考虑。最优潮流作为满 足这一目标的重要手段,近年来获得了飞速发展。
研究现状
现阶段已有的最优潮流计算方法:
• 1、非线性规划法 • 2、二次规划法 • 3、线性规划法 • 4、内点法 • 5、人工智能方法
10-8
0
2
4
6
8
10
12
迭代次数
9节点系统对偶间隙变化曲线
10-8
0
2
4
6
8
10
12
迭代次数
30节点系统对偶间隙变化曲线
三个系统的迭代次数分别为16、11、12次,迭代次数较少,计算时 间短,收敛特性好。
系统规模扩大时,迭代次数不会显著增加,说明算法对系统规模不 敏感。
初始点为非内点时,算法也能够收敛至最优解,说明算法对初始点 不敏感。
2
4
550.58 MW
A
-550.58 MW 177.08 MW
551 MW 178 Mvar
MVA
200 MW 100 Mvar
173.50 MW
20.73 MW
首先将不等式约束转化为等式约束:
g(x)l g
g(x)ug
l 0,u0
然后构造障碍函数,将含不等式约束的优化问题转化为只含等式约
束的问题:
r
r
obj. min. f (x) u log(lr ) u log(ur )
j1
j1
s.t. h(x) 0
g(x) u g
g(x) l g
构造拉格朗日函数: