对偶内点法

对偶内点法（共7篇）

对偶内点法 篇1

动态优化调度是研究在连续多个时间段内、机组运行状态已经确定的情况下,如何有效地调整各个机组在不同时段的有功功率输出,使得在满足各种静态约束和动态约束下,整个调度周期内的某些目标函数达到最优。由于动态优化调度问题的规模随着系统规模以及时段数的增加而迅速增长,因此提高算法的效率是迫切需要解决的问题[1]。

多年来,学者们通过不断探索,提出了各种求解优化调度问题的计算方法。其中,简化梯度法具有简单、容易实现等优点,但是收敛性较差,尤其是在接近最优点附近时收敛很慢[2];牛顿法和简化梯度法相比,收敛速度较快,具有二阶收敛性,但是也存在着数值不稳定、不能总是有效地确定迭代过程中起作用的约束集等缺点;人工智能算法对函数的性态要求比较低,近年来在最优潮流领域中得到了广泛的运用[3],但是存在着收敛速度比较慢的缺点。相比于以上算法,原对偶内点法具有迭代次数与系统规模关系不大、收敛性好、计算速度快等优点。文献[4]将原对偶内点法应用于非线性规模的无功优化调度问题,充分显示了内点法的收敛特性;文献[5]则采用线性约束网络流来描述电力系统动态经济调度问题,并用内点法进行求解,克服了传统方法难以精确处理不等式约束的弊端,使计算结果更精确;文献[6]则采用引入离散惩罚的非线性原对偶内点法求解动态无功优化模型,具有较快的计算速度。

随着电力行业节能减排政策的实施,在动态优化调度中同时将煤耗和污染物排放一起优化更具现实意义[7]。本文首先建立了同时考虑调度周期内煤耗量、污染物气体排放量最小的多目标动态优化调度模型,并利用模糊集理论使之转化为求取满意度最高的单目标问题;接着推导出了求解该模型的原对偶解耦内点算法;最后对算例结果进行了比较分析。

1 基于节能减排的多目标动态优化调度模型

1.1 多目标动态优化调度模型

1.1.1 目标函数

本文以常规燃煤机组的煤耗和污染物气体排放最小作为目标函数:

式(1,2)中:T为调度周期;Ng表示发电机台数;PGi(t)为第i台发电机组t时段的有功出力;ai,bi,ci为煤耗系数;αi,βi,γi为污染物排放系数。

1.1.2 静态等式约束

动态优化调度问题的静态等式约束即为各个时段的节点潮流平衡约束:

式中:省略下标t;Vi,θi为节点电压与相角;θij=θi-θj;Gij,Bij分别为系统的导纳矩阵的实部与虚部;PDi,QDi为节点i的有功负荷与无功负荷;PGi,QGi为节点的有功注入与无功注入。

1.1.3 静态不等式约束

式中不等式约束依次为发电机有功出力约束、发电机无功出力约束、节点电压约束、线路有功潮流约束。各个变量的上下限且分别用上划线和下划线表示。

1.1.4 动态约束

本文考虑发电机爬坡约束:

式中:Ramp Gi为发电机在相邻时段间能增加或减小的最大功率。

1.2 多目标动态优化调度的模糊模型

上述模型有2个相互竞争的目标函数,是多目标规划模型,有学者提出用模糊集理论来解决这个问题,并取得较好效果。文献[8]用模糊集理论将多目标函数和部分可伸缩的约束条件模糊化,使最优潮流问题在更加符合实际情况的模型上实现优化;文献[9]将模糊集理论应用于来水量不确定性分析,提出了一种水火电系统短期经济调度的新方法,较好地提高了系统的经济性;文献[10]用多目标模糊优化理论和动态规划法解决火电厂多目标负荷优化问题,得到科学、合理的负荷分配方案。因此,本文采用模糊集理论使多目标动态优化调度转化为单目标问题。在满足所有约束条件的前提下,发电煤耗和污染物气体排放总量越小越好,有上限而无下限,因此选用降半直线作为它们的隶属度函数[11],如图1、图2所示。

2个目标隶属度函数可分别用式(6,7)表示:

式中:f1(x),f2(x)分别为调度周期内的煤耗、污染物气体排放总量;c01,c02分别为单独以煤耗、污染物气体排放为优化目标时的结果;c01+δ01,c02+δ02分别为煤耗、污染物排放的最大可接受值。取λ为所有隶属函数中最小的隶属变量,可称之为满意度。

根据最大隶属度原理,可以将原来的多目标动态优化调度问题转化为满足所有约束条件的满意度最大化的问题。用以下非线性规划模型表示。

式中:xt∈Rn;ht:Rn→Rm;gt:Rn→Rr;g,g∈Rr;At∈Rq×n;gd,gd∈Rq。

上述模型包含T个时段,每个时段有n个变量、m个静态等式约束、r个静态不等式约束,各时段间有q个动态不等式约束,p个含有λ的不等式。多目标动态优化调度的模糊模型比传统动态优化调度模型增加了p个含有λ的不等式,λ不属于静态变量,因此求解该模型时更为复杂。

2 原对偶解耦内点算法

由于动态优化调度问题系统规模很大,而原对偶内点法具有迭代次数与系统规模关系不大、收敛性好、计算速度快等优点,因此本文采用原对偶内点法求解多目标动态优化调度模型。根据多目标动态优化调度模型的特点,对原对偶内点算法作出相应改进,通过对算法中修正方程的降阶与解耦,将动态优化调度的大规模线性修正方程的求解,等值解耦转换成小规模的动态变量求解和各时段静态变量求解,从而提高算法的求解效率。

2.1 引入满意度后的修正方程

将上述模型中含有λ的最后p个不等式写成根据原对偶内点法原理,首先在不等式约束中引入松弛变量,使其成为等式约束,同时在目标函数中引入对数壁垒函数,最后引入拉格朗日乘子,形成扩展问题的拉格朗日函数[12]:

式中:yt为各时段静态等式约束对应的拉格朗日乘子,各时段静态不等式约束对应的松弛变量为lt,ut>0,拉格朗日乘子为zt>0,wt>0,对数壁垒参数μt>0;动态约束对应的松弛变量为sl,su>0,拉格朗日乘子为yl,yu>0,对数壁垒参数μd>0;含有λ的不等式的松弛变量为ll,uu>0,拉格朗日乘子为zz>0,ww>0,对数壁垒参数。

为了表述方便,将各时段的优化变量以及与静态约束相关的拉格朗日乘子和松弛变量构成静态变量ρt=[xt,yt,zt,wt,lt,ut]T,每个时段的静态变量一共构成T组静态变量;定义与动态约束相关的拉格朗日乘子和松弛变量为动态变量ρd=[λ,yu,yl,su,sl,ww,zz,uu,ll]T。

原对偶内点法通过式(10)的最优性KKT条件形成一组非线性方程,然后采用牛顿法迭代求解。其中,线性化牛顿修正方程可以表示如下:

式(11)中:Wt、Lt推导过程参照文献[13],修正方程的维数高达(4r+m+n)T+4(q+p)维,需进行降阶。

2.2 修正方程的降阶

通过线性变换消去方程(11)中的Δzt,Δwt,Δlt,Δut,Δsu,Δsl,Δuu,Δll后得到:

各子矩阵WRT、LRT推导过程也可参照文献[13]。求出ΔρRT,ΔρRd后,可以通过线性变换求出其他变量。降阶后的修正方程维数为(m+n)T+2(q+p),大大减少了方程规模和计算量。

2.3 修正方程的解耦

由多目标动态优化调度的模型可知,各时段的静态约束和目标函数相互独立,相应地静态变量ρt相互独立,彼此不相关,从而使得式(11)中常数项Lt只与t时段的静态变量有关,与动态变量无关,且Δρt之间的关联系数矩阵为零矩阵;另一方面,动态约束使部分静态变量(各时段的发电机有功出力)与动态变量相互关联,因此式(11)中Δρt与Δρd之间的关联系数矩阵Et非0。由此可以得出,多目标动态优化调度问题的特殊耦合关系决定了式(11)的系数矩阵具有分块对角带边结构,降阶KKT方程(12)的系数矩阵也同样具有分块对角带边结构。针对这种特性,可以通过线性变换将式(12)等值简化成解耦方程[14]:

首先根据式(13)解出ΔρRd,然后分别求出各个时段的ΔρRt。解耦后,需要求解一个2(q+p)维方程和T个(m+n)维方程,计算量比直接求解一个(m+n)T+2(q+p)维方程要小得多。

2.4 相关参数设置

所有的对数壁垒参数都由式(15)求得,静态变量中原变量、对偶变量对应的步长根据式(16)计算,各动态变量对应步长αpd,αdd,αpλ,αdλ的计算方法同式(16)。

2.5 算法流程

算法流程如图3所示。

3 算例分析

根据上述模型和算法,本文对IEEE-30节点测试系统进行仿真计算与分析。

3.1 初始条件设置

调度周期取为一天,划分成24个时段,每个时段为1 h,图4为负荷波动系数曲线,假设系统有功功率日负荷波动曲线和无功功率日负荷波动曲线是相同的,并且各节点的负荷值在一天中以相同的负荷系数波动。各常规发电机组的煤耗系数ai,bi,ci和污染物气体排放系数αi,βi,γi如表1所示。本文以煤耗费用(元)来表征煤耗量,以排放气体的重量(t)来表征污染物气体的排放量。

以发电机有功功率上限的百分比来表示爬坡约束,本文取为5%。

单独优化煤耗的结果为10 502元,单独优化污染物气体排放时结果为3 660.2 t,取各自的最大可接受值为单独优化结果的1.4倍。

3.2 动态与静态优化调度结果对比分析

本文首先选取调度周期内煤耗量最小为单目标进行优化调度,图5比较了当爬坡约束为5%和未考虑爬坡约束时的27号发电机有功出力,即动态优化调度和静态优化调度时27号发电机有功出力变化的不同。

由图5可以看出:计及爬坡约束时,为了可以顺利到达12点时的负荷高峰,27号发电机组在第5时段作出了相应的调整,第5时段系统的负荷水平降低,但是此时煤耗较小的27号发电机就开始增加出力,从第6时段开始则全力爬坡,以保证在12点的高峰期能够承担更多的负荷;同样,为了应付22-24时段负荷的陡降,即使是在21-22时段负荷上升的情况下,27号发电机仍然减少了出力,并在23、24时段全速减小出力,保证系统能够顺利到达负荷的低点。

而如果未考虑爬坡约束,在12点负荷高峰前的多个爬坡时段和负荷陡降的22-24时段,27号发电机有功出力调节速率均超出了允许范围。由此可见,静态优化调度的结果显然是不可行的,只有动态优化调度的结果才是正确有效的。

3.3 双目标优化与单目标优化结果对比分析

图6显示了双目标优化时,各时段的煤耗量与污染物气体排放量都比单目标时的优化结果大。这是因为考虑多目标时,不仅要使得系统煤耗量最小,同时还要使污染物气体排放量最小,但这2个目标是相互矛盾的,污染物气体排放满意度的增加必然要以煤耗满意度的下降来换取,反之亦然,因此只能使得两者的综合满意度达到最高。

表3比较了双目标与单目标时整个调度周期内的总煤耗量和污染物气体总排放量。

从表3可以看到,双目标优化之后,煤耗量比单独优化煤耗时增加了256元,相对增长2.44%,但比单独优化排放时的煤耗量减少了498元,相对减少4.42%;同样,双目标优化之后,污染物气体排放量比单独优化排放时增加了89.4 t,也仅相对增长2.44%,但比单独优化煤耗时的排放量减小了461.3t,相对减小10.95%。

由此可见,虽然调度周期内的总煤耗量和排放量都没有单目标模型时的理想,但从综合效益的角度看,基于节能减排的多目标动态优化调度模型具有明显的优势,它能够更加合理的协调各个目标之间的关系。

4 结束语

本文以节能、减排作为动态优化调度的双重目标,建立了多目标动态优化调度模型,并根据模糊集理论的最大隶属度原则使之转化为求解满意度最大的单目标问题,提出了适合求解多目标动态优化问题的原对偶解耦内点算法。IEEE-30节点算例结果表明:

(1)本文通过对原对偶内点法中修正方程的降阶与解耦,大幅度提高了算法的求解效率,收敛特性较好,数值鲁棒性高。

(2)本文提出的基于节能减排的多目标动态优化调度模型能够更好的协调节约能源、减少环境污染之间的关系,并能根据机组的煤耗特性、排放特性以及爬坡约束等更加合理的安排发电机组的出力,在整个调度周期内实现各个目标函数综合满意度的最大化。

(3)本文所提的模型与算法不仅适合于双目标动态优化调度问题,还能根据实际情况灵活扩展到更多的目标函数。

基于内点法的控制分配算法及应用 篇2

关键词：控制分配,飞行控制,内点法,直接分配法,线性规划

为提高机动性和生存能力,现代高性能飞机多采用先进高效的多操纵面气动布局,这在很大程度上增加了飞控系统的控制冗余度,提高了飞机的控制能力。与传统的三种操纵面控制相比,多操纵面布局为飞机的飞行控制提供了更灵活、更可靠、更有效的实现方式。由于控制量个数大大超过了飞机的操纵输入量数目,因此如何解决多操纵面的综合分配与协调控制成为飞控系统设计中面临的首要问题[1]。

多操纵面布局飞机的出现带来了以下两个问题:

1)如何最大限度地发挥多舵面飞机的优势,将控制指令合理地分配到各操纵面上,以便提高飞机的可靠性和机动性;

2)若某操纵面在控制过程中出现饱和或发生故障,如何利用其余操纵面进行补偿,以控制飞机继续完成飞行任务。

上世纪九十年代,美国学者最早提出了控制分配技术来解决上述问题。作为一种先进的冗余控制设计方法,控制分配技术已成功地运用到了多操纵面飞控系统的设计中,并在多款机型中进行了多种控制分配方法的飞行验证[2,3,4]。

同时,控制分配技术可实现基于多目标的优化设计。针对不同的飞行任务,通过控制分配优化设计,可实现操纵面偏转能量损耗最小、雷达反射面积最小、巡航状态下阻力最小、起飞着陆状态下升力最大等不同的优化目标。

目前控制分配的算法可归为两大类:非优化分配法和优化分配法。其中,非优化控制分配方法主要是串接链法[5];基于优化的控制分配算法主要有直接分配法[6,7]、广义逆法[8,9](权值伪逆法、再分配伪逆法、多级广义逆法等)、基于线性规划[10,11](单纯形法、内点法等)和基于二次规划[12,13,14,15]的控制分配算法(序列二次规划法、定点法、有效集法、动态分配法、模型预测动态分配法等)。

线性规划中的内点法是一种多项式时间算法,具有计算效率高、收敛性好等优点[16]。本文讨论了线性规划形式的控制分配问题,采用原-对偶路径跟踪内点法[17,18]进行求解。通过仿真实验验证了该算法的有效性,并与直接分配法进行了比较。

1 控制分配问题描述

基于控制分配的飞行控制系统框图如图1所示。

其中,r为飞机的期望运动目标,设计基本控制律(设计方法通常采用动态逆法、回推法等)得到控制分配所需的伪控制量v(通常为飞机三个轴向的气动力矩或角加速度)。然后选用适当的控制分配方法,将伪控制量v合理地分配到各个作动器中去,得到每个舵机所需偏转的角度u,δ为控制舵面的偏转量,x为飞机的状态变量,y为飞机的输出。

1.1 控制分配问题的数学描述

基本的控制分配问题可描述为受约束的线性映射问题:

v∈Rn为期望的伪控制量,B∈R n×m为控制效率矩阵,行满秩(m>n),u∈R m为实际的控制量,受到方程(2)的约束。

上述问题属于一个过驱动的约束控制问题,可能有多组解,也可能只有唯一解,还可能由于约束的限制而无解。基于最优化的控制分配问题可以找到这样一个最优解:如果有多组解,找出其中最优的一个解;如果无解,找到一个使得Bu最为逼近v的解。

对于图1所示的飞行控制系统中的控制分配问题,非线性六自由度的飞机模型可用如下的方程简化描述:

作动器的动力学模型可描述为:

其中,x为飞机的状态变量,u为控制分配模块的输出,δ为作动器的输出,受(5)式的位置和速率约束。

在进行控制分配问题讨论时,通常做如下两点假设[19]:

1)由于作动器的动态响应要比飞机本身的动态响应快得多,因此在控制分配问题中常常将作动器的动力学特性忽略,即假定δ=u。

2)认为舵面偏转只产生气动力矩,只影响飞机的角加速度,忽略舵面偏转对气动力的影响。

现代飞行控制系统采用数字控制系统,因此舵机的速率约束可以转化为位置约束,因此(5)式的舵机位置和速率约束可以统一表示为舵机的位置约束:

T为采样时间。

1.2 控制分配问题的优化目标

基于优化的控制分配算法,其优化目标可分为以下三类:

(1)误差最小目标

该控制分配的优化目标是求取控制量在约束范围内的最小分配误差或加权最小误差问题,得到与期望的目标值最为接近的分配结果:

其中,vW∈R n×n为权值矩阵。

(2)控制量最小目标

仅仅满足分配误差最小并未充分发挥系统的冗余潜力,满足最小分配误差的解往往不唯一,因此又有控制量最小化目标:

(3)混合最优化目标

混合最优化问题将分配误差最小和控制量误差最小问题进行线性加权,综合考虑:

其中,h>0为加权系数,u0为期望的舵面偏转位置向量。

2 原-对偶路径跟踪内点法

Karmarkar于1984年提出了内点法用于求解线性规划问题[20],该算法不同于单纯形法,每一步迭代总是在可行解域内朝着最优解进行搜索。该算法是一种多项式时间算法,对问题的规模不敏感,具有计算量小,收敛性好的优点。它的计算复杂度为O(n3.5 L2),其中n为变量维数,L为输入长度。原-对偶路径跟踪内点法在保持解的原始可行性与对偶可行性的同时,沿着原-对偶路径找到最优解[17]。

2.1 控制分配问题转化为线性规划[10,11]

为使用线性规划求解控制分配问题,首先将控制分配问题转化为线性规划形式。考虑(10)式的混合优化目标,目标函数取l1范数,则(10)式所描述的控制分配问题可转换为如下的线性规划形式:

令,一旦获得了最优解x,则控制量为:

2.2 最优解存在的条件

考虑如下线性规划问题:

其对偶问题为:

λ为对偶变量,s为松弛变量。

由KKT(Karush-Kuhn-Tucker)条件知,若(x,λ,s)为全局最优解,当且仅当下述条件成立:

其中X=diag(x)。

引入拉格朗日乘子,则方程(13)可变为:

其中,µ>0称为障碍参数,原问题的拉格朗日乘子λ也就是对偶问题中的变量,反之亦然。

(16)式分别对x、λ求偏导,并代入e=[1 1…1]T、X=diag(x)得:

令s=µX-1e,整理上式得:

这在形式上非常接近于KKT最优解条件的方程(15),如果µ取很小的正数,则方程(18)的解几乎就是最优解。

2.3 中心路径

原-对偶可行解内点集为:

定理已经证明[16],若Ω非空,则对每个µ>0,KKT条件(15)存在唯一内点解(x(µ),λ(µ),s(µ))。将点集{(x(µ),λ(µ),s(µ))}称为原-对偶中心路径。

原-对偶路径跟踪法在求最优解的过程中,通过逐渐减小障碍参数µ的值,在原始问题和其对偶问题最优解的中心路径上移动,当µ趋于0时,原问题和对偶问题均趋于最优值。

其中0<δ<1,通常取δ=0.1。

2.4 移动方向和步长

若能求出中心路径,求取µ趋于0时的极限,即可获得最优解。然而,中心路径的解析表达式难以获得,因此,路径跟踪法并不计算中心路径,而是通过迭代大致沿着中心路径逼近最优解[16]。

任取一点(x,λ,s),其中x>0,s>0,求取一个方向(∆x,∆λ,∆s),使得新的迭代点(x+∆x,λ+∆λ,s+∆s)位于原-对偶中心路径上,即满足方程(18)式:

对上式进行整理,并忽略二阶项后,用矩阵表示为:

其中S=diag(s)。

解方程(22)得到移动方向(∆x,∆λ,∆s)。

求出移动方向后,需要确定沿此方向的步长α。为保证下一步迭代点(x+α∆x,λ+α∆λ,s+α∆s)仍然在内点集内,迭代步长α的取值应满足:

原变量与对偶变量可分别设计不同的步长αP、αD,由(23)式得到各变量的更新公式为[10]:

其中,

0

2.5 原-对偶路径跟踪内点法的算法步骤

(1)将控制分配问题转化为线性规划形式,见(11)、(12)式。

(2)给定初始点(x(1),λ(1),s(1)),其中x(1)>0,s(1)>0设定p,精度要求ε>0,正数M<∞,k=1.

(3)计算ρ=b-A x(k),σc=-A Tλ(k)-s(k),。

(4)若║ρ║1<ε,║σ║1<ε,x(k)T s(k)<ε同时成立,则停止计算,得到最优解(x(k),λ(k),s(k));若║x(k)║∞>M或║λ(k)║∞>M,则停止计算,原问题或对偶问题无界;否则进行下一步。

(5)求解移动方向和步长。

由(22)-(26)式,得(∆x(k),∆λ(k),∆s(k))及αP(k)、αD(k)。

(6)进行迭代计算:

置k=k+1,转向第(3)步。

3 仿真实验

为验证原-对偶路径跟踪内点法解决控制分配问题的有效性,以瑞典的小型单发动机战斗机模型—Admire无人机模型[19]作为对象进行控制分配仿真实验,并与Durham提出的直接分配法[6,7]进行了比较。该模型包含四组舵面,分别是鸭翼δc,右副翼δre,左副翼δle以及方向舵δr。

仿真所采用的飞行条件为0.23马赫、高度为10000英尺,30o迎角,此时的控制效率矩阵为:

舵面的位置约束为:

速率约束为:

给定的伪控制量v为该飞行点处的力矩系数,即v=(Cl,Cm,Cn)T。仿真时间为10秒,采样时间为0.02秒。

在以下的控制分配仿真结果图中,原-对偶路径跟踪内点法和直接分配法分别以ip、dir表示。图2给出了分别采用内点法和直接分配法控制分配所得到的总的命令跟踪结果。由于给定的伪控制量是一个不可达集,因此控制分配的结果与给定值存在一定误差。

两种方法的性能指标比较结果如表1所示。从表中可以看出,内点法在控制命令跟踪方面优于直接分配法;与直接分配法相比,该内点法在运算时间方面,显示出较大优势。图3给出了两种算法的时间特性比较,从图中可以清楚地看出,内点法作为一种多项式时间算法,在求解速度方面大大优于直接分配法。

图4、图5分别给出了两种控制分配方法的具体结果,即各个舵面的偏转角度和偏转角速率变化情况,相应的舵面位置约束和速率约束在图中以虚线表示。从图中可以看出,控制分配算法将伪控制指令分配到各个舵面中去,舵面的位置偏转量和偏转角速率均在舵面的约束范围内。从图4可以看出,与内点法相比,直接分配法所需的控制量更大。同时,图5中直接分配法的控制分配结果,舵面偏转存在剧烈变化,容易造成飞机飞行时的不稳定,甚至导致舵面损伤。

4 结论

为解决多操纵面布局飞机的控制指令分配问题,首先将飞行控制分配问题转化为线性规划问题,然后采用内点法中的原-对偶路径跟踪法进行求解。通过某飞机模型的仿真实验,验证了该算法的有效性,并与直接分配法进行了比较。本文飞机模型的仿真结果表明,该内点法在运算时间、分配效果方面优于直接分配法。

对偶内点法 篇3

自从20世纪60年代Carpentier提出最优潮流(OPF)问题以来,该问题受到越来越多的关注[1,2,3,4,5,6,7]。Dommel和Tinney建立了OPF问题的非线性规划模型。实际应用中,OPF有不同的表述形式,学者们也提出了不同的算法求解该问题。近年来,内点法及其改进算法由于其突出的计算性能在OPF问题中得到了广泛的应用[8,9,10,11,12,13,14]。

OPF问题中,网络潮流约束可表示为直角坐标或极坐标形式。文献[15]比较了直角坐标和极坐标2种形式的OPF模型,指出二者有相似的收敛特性和计算效率。但在OPF的内点算法中,直角坐标系的表达更简洁[16]。在潮流计算时,系统潮流方程通常表示为功率失配形式,而电流失配型潮流方程突出了电力网络本身是线性、而节点注入是非线性的特点[17],一般来说更加适合求解负荷潮流问题[18]。应用非线性内点法,文献[19]在OPF计算中采用了电流失配型潮流约束,利用电流型潮流方程一阶导数大多为常数、而二阶导数大多为0的优点,减少了编程和计算量。电流型潮流方程简化了零注入功率节点的计算,但使得其他节点的计算变得复杂。

考虑到电流和功率型潮流方程的优缺点,本文在内点法OPF问题中采用一种功率—电流混合型潮流约束模型,减少了形成雅可比矩阵和海森矩阵的计算量。多个算例的测试结果证明所提出的混合型潮流约束OPF具有更好的计算速度和收敛特性,其雅可比矩阵和海森矩阵的形成更省时、更有效;内点法的计算效率和鲁棒性得到提高,对于含零注入功率节点较多的大规模系统,其性能提升尤为明显。

1 直角坐标下的OPF问题

OPF问题是一个大规模的非凸、非线性规划问题,一般可表示为:

$\{\begin{array}{l}\underset{\mathit{x}}{{\displaystyle \min }}f(\mathit{x})\\ \text{s}.\text{t}.\mathit{h}(\mathit{x})=0\\ \underset{¯}{\mathit{g}}\le \mathit{g}(\mathit{x})\le \overline{\mathit{g}}\end{array}(1)$

式中:目标函数f(x)通常是一个具体的经济运行指标,如发电机燃料费用、系统网损等;x可以是发电机有功与无功输出、各节点电压实部和虚部、变压器分接头、移相器设置、可变并联电容器或灵活交流输电系统的控制变量等;等式约束h(x)通常表示系统潮流约束;不等式约束g(x)包括非线性不等式约束和变量的上下限约束;$\underset{¯}{\mathit{g}}$和$\overline{\mathit{g}}$分别表示其上限和下限。

1.1 目标函数f(x)

不失一般性,本文采用全系统的发电燃料费用最小为目标函数验证该功率—电流混合模型的有效性。目标函数可表示为:

$\underset{\mathit{x}}{{\displaystyle \min }}f(\mathit{x})={\displaystyle \sum _{i=1}^{n{}_{\text{G}}}(}c{}_{2i}{\rm P}{}_{\text{G}i}^2+c{}_{1i}{\rm P}{}_{\text{G}i}+c{}_{0i})(2)$

式中:c2i,c1i,c0i为第i台发电机的有功费用系数;PGi为第i台发电机的有功输出;nG为系统中可调发电机的台数。

1.2 等式约束h(x)

OPF模型中,等式约束一般为潮流方程。直角坐标下,潮流约束方程可表示为功率失配和电流失配2种形式。

功率失配形式如下:

$\{\begin{array}{l}{\rm P}{}_{\text{G}\text{D}i}-e{}_i{\displaystyle \sum _{j\in i}(}G{}_{ij}e{}_j-B{}_{ij}f{}_j)-f{}_i{\displaystyle \sum _{j\in i}(}G{}_{ij}f{}_j+\\ B{}_{ij}e{}_j)=0\\ Q{}_{\text{G}\text{D}i}-f{}_i{\displaystyle \sum _{j\in i}(}G{}_{ij}e{}_j-B{}_{ij}f{}_j)+e{}_i{\displaystyle \sum _{j\in i}(}G{}_{ij}f{}_j+\\ B{}_{ij}e{}_j)=0\end{array}(3)$

式中:i=1,2,…,nb;nb为网络节点总数;ei和fi分别为节点i电压的实部和虚部;Gij和Bij分别为系统节点导纳矩阵的实部和虚部。

对于接有发电机或负荷的系统节点,PGDi和QGDi可分别表示为:

$\{\begin{array}{l}{\rm P}{}_{\text{G}\text{D}i}={\displaystyle \sum _{j\in i}{\rm P}}{}_{\text{G}j}-{\displaystyle \sum _{k\in i}{\rm P}}{}_{\text{D}k}\\ Q{}_{\text{G}\text{D}i}={\displaystyle \sum _{j\in i}Q}{}_{\text{G}j}-{\displaystyle \sum _{k\in i}Q}{}_{\text{D}k}\end{array}(4)$

式中:PGj和QGj分别为有功和无功输入;PDk和QDk分别为有功和无功负荷。

对于没有连接任何发电机和负荷的节点,PGDi=QGDi=0。

直角坐标下,电流失配型潮流方程可表示为:

$\{\begin{array}{l}{\rm I}{}_{xi}-{\displaystyle \sum _{j\in i}(}G{}_{ij}e{}_j-B{}_{ij}f{}_j)=0\\ {\rm I}{}_{yi}-{\displaystyle \sum _{j\in i}(}G{}_{ij}f{}_j+B{}_{ij}e{}_j)=0\end{array}(5)$

式中:i=1,2,…,nb;Ixi和Iyi分别为节点i注入电流的实部和虚部,对于接有发电机或负荷的系统节点,Ixi和Iyi可表示为:

$\{\begin{array}{l}{\rm I}{}_{xi}=\frac{{\rm P}{}_{\text{G}\text{D}i}e{}_i+Q{}_{\text{G}\text{D}i}f{}_i}{e{}_i^2+f{}_i^2}\\ {\rm I}{}_{yi}=\frac{{\rm P}{}_{\text{G}\text{D}i}f{}_i-Q{}_{\text{G}\text{D}i}e{}_i}{e{}_i^2+f{}_i^2}\end{array}(6)$

对于没有连接任何发电机和负荷的节点,Ixi=Iyi=0。

1.3 不等式约束g(x)

OPF模型中,需考虑优化运行点的静态安全约束。本文考虑如下的不等式约束:

$\{\begin{array}{l}{\rm P}{}_{\text{G}i\text{m}\text{i}\text{n}}\le {\rm P}{}_{\text{G}i}\le {\rm P}{}_{\text{G}i\text{m}\text{a}\text{x}}1\le i\le n{}_{\text{G}}\\ Q{}_{\text{G}i\text{m}\text{i}\text{n}}\le Q{}_{\text{G}i}\le Q{}_{\text{G}i\text{m}\text{a}\text{x}}1\le i\le n{}_{\text{G}}\\ V{}_{i\text{m}\text{i}\text{n}}^2\le e{}_i^2+f{}_i^2\le V{}_{i\text{m}\text{a}\text{x}}^{2}1\le i\le n{}_{\text{b}}\\ t{}_{k\text{m}\text{i}\text{n}}\le t{}_k\le t{}_{k\text{m}\text{a}\text{x}}1\le k\le n{}_{\text{t}}\\ S{}_k^2\le S{}_{k\text{m}\text{a}\text{x}}^{2}1\le k\le n{}_{\text{l}}\end{array}(7)$

式中:S2k=P2k+Q2k;nt为可调变压器的台数;nl为电力网络中功率受限的支路条数。

本文考虑的OPF模型中,x=[PG,QG,t,e,f]T,其中,PG,QG,t为控制变量,e和f为状态变量。

2 潮流方程的混合模型

考虑到功率失配和电流失配型潮流方程各自的优缺点,本文提出了一种功率—电流混合型潮流约束模型。在直角坐标下,该混合模型将所有节点分成2类:非零注入功率节点(用集合I1表示)和零注入功率节点(用集合I2表示)。考虑到实际电力系统的结构和特点,本文定义零注入功率节点如下:

1)没有连接发电机、负荷及补偿装置的联络节点是零注入功率节点;

2)对于装设有补偿装置的节点,由于补偿装置一般为电容和电感,一般视为阻抗元件,该阻抗可包含在节点导纳矩阵中,这样的节点可视为零注入功率节点;

3)对于部分可视为恒定阻抗的负荷,则该负荷可以用一恒定阻抗模拟,从而包含于节点导纳矩阵中,则仅接有这类负荷的节点也可以视为零注入功率节点。

在图1所示的简单电力系统中,按照上述定义,节点2为零注入功率节点,实际电力系统中,可视为零注入功率的节点数量非常多。仅那些接有发电机(节点1)或不能等效为恒定阻抗负荷(节点3)的节点被视为非零注入功率节点,实际电力系统中这类节点较少。实际电力系统的这种特点有利于应用本文提出的混合模型,充分挖掘混合模型的优势。

该混合模型中,对于零注入功率节点和非零注入功率节点,其潮流方程分别表示如下。

对∀i∈I1,

$\{\begin{array}{l}{\rm P}{}_{\text{G}\text{D}i}=e{}_i{\displaystyle \sum _{j\in i}(}G{}_{ij}e{}_j-B{}_{ij}f{}_j)+\\ f{}_i{\displaystyle \sum _{j\in i}(}G{}_{ij}f{}_j+B{}_{ij}e{}_j)\\ Q{}_{\text{G}\text{D}i}=f{}_i{\displaystyle \sum _{j\in i}(}G{}_{ij}e{}_j-B{}_{ij}f{}_j)-\\ e{}_i{\displaystyle \sum _{j\in i}(}G{}_{ij}f{}_j+B{}_{ij}e{}_j)\end{array}(8)$

对∀i∈I2,

$\{\begin{array}{l}{\displaystyle \sum _{j\in i}(}G{}_{ij}e{}_j-B{}_{ij}f{}_j)=0\\ {\displaystyle \sum _{j\in i}(}B{}_{ij}e{}_j+G{}_{ij}f{}_j)=0\end{array}(9)$

在基于混合型潮流约束的OPF模型中,潮流约束方程h(x)=0可对应地分解为2部分:h1(x)=0,由式(8)表示;h2(x)=0,由式(9)表示,即

$\mathit{h}(\mathit{x})=\left[\begin{array}{l}\mathit{h}{}_1(\mathit{x})\\ \mathit{h}{}_2(\mathit{x})\end{array}\right](10)$

注意到混合模型(10)具有以下特点:对于节点i∈I2,其一阶导数为常数,对应的雅可比矩阵元素不变[18];对于节点i∈I2,其二阶导数为0,对应的海森矩阵元素为0。因此,可以肯定集合I2中的所有节点对拉格朗日函数的海森矩阵没有影响,且这些节点对应的雅可比矩阵元素为恒定值,可在初始化阶段得到。

3 内点法中功率—电流混合潮流约束处理

应用内点法求解非线性规划问题(式(1))时,首先需构造如下形式的拉格朗日函数[20]:

$\begin{array}{l}L{}_{\text{g}}=f(\mathit{x})-\mathit{y}{}^{\text{Τ}}\mathit{h}(\mathit{x})-\mathit{z}{}^{\text{Τ}}(\mathit{g}(\mathit{x})-\mathit{l}-\underset{¯}{\mathit{g}})-\\ \mathit{w}{}^{\text{Τ}}(\mathit{g}(\mathit{x})+\mathit{u}-\overline{\mathit{g}})-\mu {\displaystyle \sum _{j=1}^r\ln }l{}_j-\mu {\displaystyle \sum _{j=1}^r\ln }u{}_j(11)\end{array}$

式中:l和u为松弛向量;μ为障碍参数;y,z,w为对应等式约束和不等式约束的拉格朗日乘子;r为不等式约束维数。

上述拉格朗日函数的KKT一阶最优条件为:

$\{\begin{array}{l}\mathit{L}{}_{\mathit{x}}=\nabla {}_{\mathit{x}}f(\mathit{x})-\nabla {}_{\mathit{x}}\mathit{h}(\mathit{x})\mathit{y}-\nabla {}_{\mathit{x}}\mathit{g}(\mathit{x})(\mathit{z}+\mathit{w})=0\\ \mathit{L}{}_{\mathit{y}}=\mathit{h}(\mathit{x})=0\\ \mathit{L}{}_{\mathit{z}}=\mathit{g}(\mathit{x})-\mathit{l}-\underset{¯}{\mathit{g}}=0\\ \mathit{L}{}_{\mathit{w}}=\mathit{g}(\mathit{x})+\mathit{u}-\overline{\mathit{g}}=0\\ \mathit{L}{}_{\mathit{l}}=\mathit{L}\mathit{{\rm Z}}\mathit{E}-\mu \mathit{E}=0\\ \mathit{L}{}_{\mathit{u}}=\mathit{U}\mathit{W}\mathit{E}+\mu \mathit{E}=0\end{array}(12)$

式中:L=diag(l1,l2,…,lr);U=diag(u1,u2,…,ur);Z=diag(z1,z2,…,zr);W=diag(w1,w2,…,wr);E=[1,1,…,1]T。

采用牛顿法求解上述非线性方程,可得到如下3个解耦的线性方程组:

$\left[\begin{array}{cc}{\mathit{{\rm H}}}^{\prime }& \nabla {}_{\mathit{x}}\mathit{h}(\mathit{x})\\ \nabla {}_{\mathit{x}}^{\text{Τ}}\mathit{h}(\mathit{x})& 0\end{array}\right]\left[\begin{array}{c}\text{Δ}\mathit{x}\\ \text{Δ}\mathit{y}\end{array}\right]=\left[\begin{array}{c}\mathit{L}{}_{{\mathit{x}}^{\prime }}\\ -\mathit{L}{}_{\mathit{y}}\end{array}\right](13)$

$\left[\begin{array}{cc}\mathit{L}& \mathit{{\rm Z}}\\ 0& \mathit{{\rm I}}\end{array}\right]\left[\begin{array}{c}\text{Δ}\mathit{z}\\ \text{Δ}\mathit{l}\end{array}\right]=\left[\begin{array}{c}-\mathit{L}{}_{{\mathit{l}}^{\prime }}\\ \mathit{L}{}_{\mathit{z}}+\nabla {}_{\mathit{x}}^{\text{Τ}}\mathit{g}(\mathit{x})\text{Δ}\mathit{x}\end{array}\right](14)$

$\left[\begin{array}{cc}\mathit{U}& \mathit{W}\\ 0& \mathit{{\rm I}}\end{array}\right]\left[\begin{array}{c}\text{Δ}\mathit{w}\\ \text{Δ}\mathit{u}\end{array}\right]=\left[\begin{array}{c}-\mathit{L}{}_{{\mathit{u}}^{\prime }}\\ -\mathit{L}{}_{\mathit{w}}-\nabla {}_{\mathit{x}}^{\text{Τ}}\mathit{g}(\mathit{x})\text{Δ}\mathit{x}\end{array}\right](15)$

式中:

Ll′=LZE-μE+diag(Δl)Δz

Lu′=UWE+μE+diag(Δu)Δw

引入预测—校正内点法,采用如下步骤求解该非线性规划问题:

1)初始化:置迭代次数k=0,给定控制变量和状态变量、松弛变量、拉格朗日乘子的初值x0,l0,u0,y0,z0,w0。

2)计算目标函数和约束的雅可比矩阵和海森矩阵,形成式(13)～式(15)所示的线性系统。

3)对式(13)的系数矩阵进行LU分解。

4)令μ=0,diag(Δl)Δz=0,diag(Δu)Δw=0,求解式(13),得到预估值Δxaf,Δyaf,从而通过回代求解式(14)和式(15),得到Δlaf,Δuaf,Δzaf,Δwaf。

5)根据下式计算仿射步长αafp和αafd:

$\{\begin{array}{l}\alpha {}_{\text{a}\text{f}\text{p}}=\min \left(0.9995\underset{i}{{\displaystyle \min }}\left(\frac{-l{}_i}{\text{Δ}l{}_{\text{a}\text{f}i}},\frac{-u{}_i}{\text{Δ}u{}_{\text{a}\text{f}i}}\right),1\right)\\ \text{Δ}l{}_{\text{a}\text{f}i}<0,\text{Δ}u{}_{\text{a}\text{f}i}<0\\ \alpha {}_{\text{a}\text{f}\text{d}}=\min \left(0.9995\underset{i}{{\displaystyle \min }}\left(\frac{-z{}_i}{\text{Δ}z{}_{\text{a}\text{f}i}},\frac{-w{}_i}{\text{Δ}w{}_{\text{a}\text{f}i}}\right),1\right)\\ \text{Δ}z{}_{\text{a}\text{f}i}<0,\text{Δ}w{}_{\text{a}\text{f}i}>0\end{array}(16)$

更新中心参数σ和障碍参数μaf:

$\{\begin{array}{l}\sigma =\min \left(\left(\frac{\rho {}_{\text{a}\text{f}}}{\rho }\right){}^3,0.1\right)\\ \mu {}_{\text{a}\text{f}}=\sigma \frac{\rho }{2r}\end{array}(17)$

式中:ρ=lTz-uTw;ρaf=(l+αafpΔlaf)T(z+αafdΔzaf)-(u+αafpΔuaf)T(w+αafdΔwaf)。

6)令μ=μaf,diag(Δl)Δz=diag(Δlaf)Δzaf,diag(Δu)Δw=diag(Δuaf)Δwaf, 采用步骤3所得LU矩阵求解式(13),得到总的预测—校正修正量Δx,Δy;通过回代求解式(14)和式(15)求得Δl,Δu,Δz,Δw;按式(16)更新步长αp和αd,其中下标为af的变量分别用本步骤所得的变量代替;按下式更新原对偶变量:

$\{\begin{array}{l}\mathit{x}=\mathit{x}+\alpha {}_{\text{p}}\text{Δ}\mathit{x}\\ \mathit{l}=\mathit{l}+\alpha {}_{\text{p}}\text{Δ}\mathit{l}\\ \mathit{u}=\mathit{u}+\alpha {}_{\text{p}}\text{Δ}\mathit{u}\\ \mathit{y}=\mathit{y}+\alpha {}_{\text{d}}\text{Δ}\mathit{y}\\ \mathit{z}=\mathit{z}+\alpha {}_{\text{d}}\text{Δ}\mathit{z}\\ \mathit{w}=\mathit{w}+\alpha {}_{\text{d}}\text{Δ}\mathit{w}\end{array}(18)$

7)计算对偶间隙ρ。若对偶间隙和最大潮流偏差的绝对值小于给定精度或达到最大迭代次数,则停止计算,否则转步骤2。

第2节已经指出,直角坐标下的功率—电流混合型潮流约束中,所有节点可归入非零注入功率节点集I1和零注入功率节点集I2。对于节点i∈I1,采用功率失配型潮流方程;对于节点i∈I2,采用电流失配型潮流方程。

节点导纳矩阵G和B也按上述分类分成2部分:

$\{\begin{array}{l}\mathit{G}=\left[\begin{array}{l}\mathit{G}{}_1\\ \mathit{G}{}_2\end{array}\right]\\ \mathit{B}=\left[\begin{array}{l}\mathit{B}{}_1\\ \mathit{B}{}_2\end{array}\right]\end{array}(19)$

内点法求解过程中每次迭代需要求解xh(x),xh(x)y和2xh(x)y。在混合模型中,xh(x)也可相应地分成2部分:

$\nabla {}_{\mathit{x}}\mathit{h}(\mathit{x})=\left[\begin{array}{l}\nabla {}_{\mathit{x}}\mathit{h}{}_1(\mathit{x})\\ \nabla {}_{\mathit{x}}\mathit{h}{}_2(\mathit{x})\end{array}\right](20)$

若变压器分接头不可调,则x=[PG,QG,e,f]T,其中仅PG和QG为控制变量,e和f为状态变量。I2对应的雅可比矩阵和海森矩阵为:

$\{\begin{array}{l}\nabla {}_{\mathit{x}}\mathit{h}{}_2(\mathit{x})=\left[\begin{array}{cccc}0& 0& \mathit{G}{}_2& -\mathit{B}{}_2\\ 0& 0& \mathit{B}{}_2& \mathit{G}{}_2\end{array}\right]\\ \nabla {}_{\mathit{x}}^2\mathit{h}{}_2(\mathit{x})=0\\ \nabla {}_{\mathit{x}}^2\mathit{h}{}_2(\mathit{x})\mathit{y}=0\end{array}(21)$

这意味着h2(x)的一阶导数为常数,而二阶导数为0。因此,h2(x)对拉格朗日函数Lg的海森矩阵元素没有影响,从而可减少形成雅可比矩阵和海森矩阵的计算量,即

2xh(x)y=2x(y${}_1^{\text{Τ}}$h1(x)+y${}_2^{\text{Τ}}$h2(x))=

2x (y${}_1^{\text{Τ}}$h1(x)) (22)

其中,等式约束的拉格朗日乘子也对应地分成2部分,即y=[y1,y2]T。

若变压器分接头可调,则x=[PG,QG,t,e,f]T,此时PG,QG,t均为控制变量,且h2(x)中t对应的节点对拉格朗日函数的海森矩阵中少数元素有影响,但其雅可比矩阵和海森矩阵的迭代更新计算量仍然少于全部采用功率失配模型的内点法OPF模型。与纯功率失配型潮流约束相比,本文所提出的混合型潮流约束能够减少零注入功率节点的计算量以及内点算法中修正方程的系数矩阵非零元个数,从而减少线性系统的计算时间。

如果电力系统含较多的零注入功率节点,则形成雅可比矩阵和海森矩阵的计算量将明显减少,且线性方程组(13)的求解速度也得到提升。该混合型潮流约束能够改善内点法OPF问题的收敛特性,对于含大量零注入功率节点的电力系统OPF问题效果尤为明显。

4 算例分析

本文采用了6个较大规模算例测试所提出的混合模型在OPF问题中的效果。表1给出了这6个算例的参数,其中,2 052节点系统数据取自日本某实际电力系统,其余测试系统的参数均取自电力系统计算软件MATPOWER,所有补偿装置均视为恒电抗,且认为变压器变比连续可调。为求解大规模问题,测试中充分利用了稀疏矩阵技术。对于所有算例,最大迭代次数均取为100,迭代精度设置为对偶间隙和最大潮流方程偏差均小于10-6,软件测试环境为Microsoft Visual C++2005,硬件条件为Intel 1.73 GHz主频,1 GB内存。

为比较混合型、功率型和电流型潮流方程的收敛特性和计算效率,本文对于各变量均采用平启动策略,即各控制变量、状态变量、松弛变量和拉格朗日乘子按下式取初值:

$\begin{array}{l}\{\begin{array}{l}{\rm P}{}_{\text{G}i}^0=\frac{{\rm P}{}_{\text{G}i\text{m}\text{i}\text{n}}+{\rm P}{}_{\text{G}i\text{m}\text{a}\text{x}}}{2}\\ Q{}_{\text{G}i}^0=\frac{Q{}_{\text{G}i\text{m}\text{i}\text{n}}+Q{}_{\text{G}i\text{m}\text{a}\text{x}}}{2}\\ e{}_i^0=1.0\\ f{}_i^0=0.0\\ t{}_k^0=\frac{t{}_{k\text{m}\text{i}\text{n}}+t{}_{k\text{m}\text{a}\text{x}}}{2}\end{array}(23)\\ \{\begin{array}{l}\mathit{l}{}^0=\gamma (\overline{\mathit{g}}-\underset{¯}{\mathit{g}})\\ \mathit{u}{}^0=(1-\gamma )(\overline{\mathit{g}}-\underset{¯}{\mathit{g}})\\ z{}_i^0=\frac{\mu {}_0}{l{}_{0i}}\\ w{}_i^0=\frac{-\mu {}_0}{u{}_{0i}}\\ \mathit{y}{}^0=1.0\end{array}(24)\end{array}$

式中:γ=0.64,μ0=nb。

表2给出了所有算例给定目标函数下的迭代次数Niter和计算时间Tcpu,其中计算时间指不考虑原始数据的输入输出、形成导纳矩阵等,仅考虑迭代求解过程的耗时,包括每次迭代更新雅可比矩阵和海森矩阵以及求解KKT系统的时间。

为比较混合模型的计算效率,本文将测试算例分为含大量零注入功率节点的算例A和含少量零注入功率节点的算例B。从表1可知,2 052节点测试系统的零注入功率节点比例高达76.95%,可归为A类;其他算例中零注入功率节点的比例均小于35%,可归为B类。

2类测试算例在不同潮流约束形式下的内点算法性能比较结果显示,对于A类算例,在平启动方式下仅混合型潮流约束模型在给定的目标函数下收敛,采用纯功率型和纯电流型潮流约束模型的内点法OPF均不收敛。这说明对于含大量零注入功率节点的系统,采用混合型潮流约束方程具有更好的收敛特性和计算效率。

对于B类算例,在相同的初始条件下,功率型潮流约束内点OPF总迭代次数为104次,总计算时间为9.48 s;电流型潮流约束内点OPF总迭代次数为105次,与功率型相当,但计算时间却减少至8.65 s,这体现了电流型潮流方程的优点;而结合了功率型和电流型潮流约束优点的混合型潮流约束内点OPF总迭代次数仅为96次,总计算时间仅为7.80 s。这表明,对于含较少零注入功率节点的系统,该混合型潮流约束与功率型或电流型潮流约束相比具有更好的收敛特性和计算效率,只是效果没有含大量零注入功率节点的系统那么明显。

5 结语

本文提出了一种直角坐标下的功率—电流混合型潮流约束模型,并采用内点法求解OPF问题。该混合形式将所有节点分为零注入功率节点和非零注入功率节点2类。对于零注入节点,采用电流失配型潮流约束;对于非零注入节点,采用功率失配型潮流约束。该模型结合了功率型和电流型潮流约束二者的优点。对于最多含2 746个节点的多个不同规模的算例测试证明,该混合形式比单纯采用功率型或电流型潮流约束具有更好的收敛特性和计算效率,尤其适合求解含大量零注入功率节点系统的内点法OPF问题。

摘要：提出了一种基于直角坐标下功率—电流混合型潮流约束的最优潮流模型。该模型对系统中的非零注入功率节点采用功率失配型潮流约束,而对零注入功率节点采用电流失配型潮流约束。这种混合模型结合了功率和电流型潮流方程的优点:对于零注入功率节点,该模型具有电流型潮流方程一阶导数为常数、二阶导数为0的特点,从而使雅可比矩阵和海森矩阵更容易计算;对于非零注入功率节点,该模型也比电流型潮流方程更好处理。该模型特别适合非线性预测—校正内点法的最优潮流,多个大规模算例测试证明该模型收敛性更好,计算效率更高,尤其适合求解含较多零注入功率节点的大规模电力系统最优潮流问题。

对偶内点法 篇4

当今, 分布式电源在配电网得到广泛的运用, 快速的发展。DG有利于减少用户的电能花费, 缓解电网的拥堵, 在负荷集中点安装环保能源, 可以提高电压稳定性, 减小网络损耗, 缓解储备容量［1］。DGs一般指发电量在1 k W到50 MW之间, 安装在负荷集中区的发电电源［2］。

近年来, 国内外的大量的学者在这方面做了大量的研究。Sudipta Ghosh等采用牛顿拉夫逊求解网损和花费最小的DGs的位置, 获得最大的经济效益［2］。Luis F等采用多期交流优化潮流求解能耗最小的确定DG位置［3］。Isrsfil Hussain和Anjan采用DE方法以网损最小进行选址定容［4］。M. F. Alhajri等采用FSQP方法以网损最小进行选容［5］。但是对于多目标选址定容, 依然不能综合考虑优化。

内点法已被证明是解决非线性规划的一种强有力工具［6］, 表现出极好的收敛性和较高的精度［7-10］, 在电力系统领域上得到了广泛的应用［11-13］, 但在解决有离散问题上存在不足。在解决多目标问题上, 基于Pareto最优意义的协调各目标函数之间的关系［14］; 采用模糊理论适合描述不确定性及处理不同量纲及互相矛盾的多目标优化问题［15］, 把多目标函数转换成单目标函数, 通过模糊选择控制［16］实行。这些处理方法和人工智能算法具有很好的结合性, 对于内点法不能很好的实现。基于此, 论文从另一个角度考虑分布式电源的选址定容, 充分利用内点法收敛性好, 精度高等优势。

1 多目标模型

1. 1 目标函数

1) 以有功网损最小:

2) 节点电压水平:

式中, ︱V*i︱=1, ︱Vi︱为节点i的指定电压幅值, Δ︱Vimax︱表示节点i节点允许的最大电压偏差, NPQ为系统的负荷节点数。

把Ploss和 ΔV采用式 ( 3) 处理, μ∈ (0, 1) , 其反应目标子函数的水平。定义如下:

式中: i = 1, 2; μ1, μ2对应于系统网络损耗和节点电压水平的子目标函数。

式中: w协调因子。

1. 2 约束条件

1. 2. 1 功率方程

1.2.2不等式约束

线性不等式约束:

非线性约束条件:

nbr为支路数; nDG为DG台数; ng为发电机台数。

2 模型的分解

把DG的选址定容的数学模型写成标准形式:

式中: xc为连续变量向量, xc= [δ v PgQgSDG]T, xd为离散变量向量, DG安装的PQ节点号。

3 内点法

连续部分采用内点法如下, 即把式 ( 9) 引入松弛变量Zm转化为:

根据Karush - Kuhn - Tucker最优一阶必要条件得到:

对最优化条件式 ( 11) 采用牛顿法求解得到式 ( 12) 。

牛顿跌代更新计算可以根据以下3步:

1) 根据式 (14) 计算Δx和Δλ;

2) 根据式 (13) 计算ΔZ;

3) 根据式 (12) 计算Δμ。

αp、αd分别为原变量和对偶变量步长, 表达如式 ( 15) 和式 ( 16) :

变量更新如式 ( 17) :

式中: ξ = 0. 999 95, 在牛顿迭代中, 为满足KKT则有扰动参数 γ 趋于0。。

4 算例分析

论文采用Matpower4. 1 中的IEEE30 节点作为测试模型 ( 改进IEEE30 bus) 。以节点1 作为平衡节点, 其基准电压为100 MV。

4. 1协调因子w优化

对于case 30bus以无DG网损为目标函数 ΔVmax=30. 556, 安装2 台DG以节点电压 ΔVmin= 29. 205。以无DG节点电压为目标函数Ploss，max= 0. 025, 2 台DG损耗为目标函数Ploss，min= 0. 015 以取case 30 busΔV区间[29. 2056, 30. 556 3 ], case30 bus Ploss区间[0. 014 8, 0. 024 9]。建立多目标函数, w和Ploss, ΔV的关系如图1 所示。

为达到子目标函数优化效果 ( 主要考虑网损最小) , 选取协调因子w = 0. 4。

4. 2 DG容量和最优位置

如图2 所示, 最优位置是8 节点, P = 0. 236, Q =0. 047 2, S = 0. 242 9, 在IEEE30 节点模型中, 虽然节点7, 8不是负荷集中区, 却是负荷的最严重区, 所以DG安装8节点合理性。故安装2 台DG, 安装位置为8, 11 节点, S8= 0. 242 9, S11= 0. 156 9。

4. 3 优化的效果

如图3 所示, 对于安装1 台DG建立的多目标函数比没有安装DG以网损为目标函数的电压稳定性总体上, 有大幅度提高; 同样, 相比较没有安装DG以节点电压为目标函数, 更趋于平稳更靠近1 点附近。空白柱形表示安装有DG以多目标函数; 斜线柱形表示没有安装DG以节点电压水平为目标函数, 交叉线柱形表示没有安装DG以网损为目标函数。

采用上述方法安装1 台DG时, Ploss= 0. 016 6 比较没有安装DG以网络有功损耗为目标函数Ploss= 0. 022 7, 网损减少了36. 75% , 安装1 台DG以网络有功损耗为目标函数Ploss= 0. 015 7, ΔVloss= 30. 471 6, 网损增加了5. 42% ;ΔVloss比相同条件下以节点电压水平为目标函数求得的ΔVu= 29. 070 9, 增加了4. 8% , 故满足优化要求。 ( 一次连续迭代次数约31 次, 时间23. 462 s。)

5 总结

对偶内点法 篇5

随着电力系统的不断发展，最优潮流问题在电力系统运行和控制领域有着更为重要的意义。法国学者Carpentier[1]于20世纪60年代最先提出了严格的电力系统最优潮流数学模型，之后国内外众多专家、学者对其进行了深入的研究，主要包括非线性规划算法[2]、二次规划算法[3]、线性规划算法[4]、梯度及牛顿类算法[5]、内点算法[6,7]和人工智能算法[8,9,10]等各种优化算法。其中，美籍印度学者Karmarkar[7]提出的内点法因其具有多项式时间复杂特性，相对于其他优化算法在求解大规模优化问题上有着明显的优势，成为目前求解最优潮流问题最为广泛的算法。

电力系统最优潮流问题作为复杂的非线性规划问题，不仅要求算法具有良好的收敛性，还要求算法具有快速计算能力，以满足电力系统在线计算的要求。虽然改进后的内点算法[11,12,13]在求解大规模电网最优潮流时具有鲁棒性好、收敛性强等优点，但是在非线性约束的处理上仍然比较耗时，难以满足在线计算的要求。直流潮流（direct current power flow,DCPF)[14,15]的出现对于电力系统静态安全分析和电力市场有着重要的意义。然而，传统DCPF模型对实际系统进行了大量简化和近似，完全忽略了电压和无功功率的影响，所以计算精度并不理想。基于此，文献[15]通过建立线路有功功率和无功功率之间的数学关系式，提出了一种类DCPF算法，该算法提高了DCPF的计算精度。随着DCPF的发展，直流最优潮流（direct current optimal power flow,DCOPF)[16,17]也逐渐受到了人们的关注。与交流最优潮流（alternating current optimal power flow,ACOPF）相比，DCOPF通过对传统最优潮流问题中的非线性约束进行一定的简化和近似，最后转变为仅含线性约束的最优潮流问题，可以大幅提高计算速度，但是计算精度不高。

为了提高DCOPF的计算精度，本文提出了拟直流最优潮流（pseudo direct current optimal power flow,PDCOPF）模型，该模型利用有功功率和无功功率之间的数学联系，对有功功率平衡等式约束进行修正，从而计及了无功功率的影响。此外，通过采用简化预测—校正内点法（simplified predictorcorrector interior point method,SPCIPM）求解最优潮流问题，保证了算法的计算效率。

1 最优潮流的数学模型

一般非线性优化问题的数学模型可以描述如下：

式中：f(x）为目标函数；h(x）为等式约束；g(x）为不等式约束；gmax和gmin分别为不等式约束的上、下限。

最优潮流问题的目标函数种类很多，本文采用系统运行成本最小，即火电机组燃料费用最小为目标函数，表达如下：

式中：x=[Pg,Qg,V，θ]T，其中Pg和Qg分别为发电机的有功、无功出力，V和θ分别为节点电压幅值、相角；Pgk为发电机k的有功出力；ak,bk,ck为发电机k的耗量特性参数；ng为发电机数。

最优潮流问题中的等式约束主要为电网功率平衡方程，如下所示：

式中：Psi和Qsi分别为节点i的有功、无功注入功率；Vi和Vj分别为节点i、节点j的电压幅值；θij为节点i和节点j之间的电压相角差；Gij和Bij分别为节点导纳矩阵第i行第j列元素的实部、虚部；Pij和Qij分别为线路ij的有功、无功功率；n为节点数。

对于线路ij，在忽略其并联导纳的情况下，线路功率方程可以描述为：

式中：rij和xij分别为线路ij的电阻、电抗。

最优潮流中的不等式约束主要为变量约束或变量函数约束，如下所示：

式中：Qgk为发电机k的无功出力；Pgk,max,Pgk,min,Qgk,max,Qgk,min分别为发电机k有功出力的上、下限和无功出力的上、下限；Qi为节点i的电压相角；Vi,max,Vi,min，θi,max，θi,min分别为节点i的电压幅值上、下限和电压相角上、下限。

2 基于拟直流模型的SPCIPM

2.1基于拟直流模型的最优潮流

在正常运行的实际电力系统中，各节点电压一般稳定在额定电压附近，线路两端的电压相角差很小，并且对于超高压电力网络，线路电阻远小于线路电抗[18]。因此可做如下简化假设：Vi=Vj=1,sinθij=θij,cosθij=1,rij=0，从而有Qij=0，式（5）可以化简为：

结合式（11），式（3）写成矩阵形式：

式中：ΔP,Ps，θ均为n维列向量；B0为以1/xij为线路导纳建立起来的n×n阶节点导纳矩阵，矩阵中各元素如式（13）所示。

综上，在最优潮流的标准直流模型中，目标函数为式（2），其中，x=[Pg，θ]T；等式约束为式（12），不等式约束为式（7）和式（10）。由于标准直流模型将非线性约束转化为线性约束，可以提高计算效率，但是却忽略了电压和无功功率的影响，因此计算结果的精度不高。

基于此，首先将式（5）和式（6）联立，可得：

对式（14）做如下简化近似：Vi=Vj=1,sinθij=θij。则其可以进一步写为：

式中：αij为线路ij的电阻与电抗的比值。

结合式（15），式（3）可以改写为：

式中：Pcor为无功功率对有功功率平衡方程的修正量，其中的元素

综上，最优潮流的拟直流模型只要在标准直流模型的基础上将等式约束由式（12）变为式（16）即可。在拟直流模型中，由于引入无功功率对有功功率平衡方程进行修正，并且考虑了线路电阻，数学模型更接近交流模型，比标准直流模型有着更高的精度。

2.2 SPCIPM

文献[19]详细介绍了采用原—对偶内点法求解最优潮流的过程，由于文中对不等式上下限约束的独立处理，在采用拉格朗日乘子法求解时，需要引入多余的松弛变量和对应的拉格朗日乘子，求解过程变得更为复杂和繁琐，同时也影响了算法的效率。

基于此，本文提出了一种简化内点算法，通过对不等式约束的简化处理，形成只含上限约束的广义不等式约束，减少了松弛变量和对应拉格朗日乘子的引入，大大简化程序的同时提高了算法的效率。简化原—对偶内点法的具体求解过程[19]详见附录A。

为了进一步提高计算效率，本文采用收敛性更好的SPCIPM求解PDCOPF。该方法通过预测步和校正步来计算考虑高阶项影响的牛顿方向，从而可以获得更大的迭代步长，使牛顿方向接近中心收敛路径，进一步提高了算法的收敛性。其中主要的区别在于SPCIPM在附录A式（A7）泰勒展开过程中考虑了二阶项的影响，即将附录A式（A14）改写为：

再结合附录A式（A7），有

具体的预测和校正过程详见下文。

3 PDCOPF的计算流程

基于SPCIPM的PDCOPF的求解过程可以归纳如下[11,12]。

步骤1：初始数据设置和常量计算（导纳矩阵等），包括变量x,u和拉格朗日乘子y,w，并置迭代次数K=0，最大迭代次数Kmax=50，中心参数σ=0。

步骤2：计算互补间隙Cgap，若满足Cgap<ε（收敛精度ε一般取10-8～10-6），则结束程序并输出最优解；否则继续。

步骤3：根据每次迭代中的发电机有功出力，进行传统交流潮流计算，求得对应的线路无功潮流，以修正有功功率平衡约束方程。

步骤4：计算雅可比矩阵"xh(x),"xg～(x）和各常数项Lx′，Ly,Lw,Lu，μ。

步骤5：预测步。根据附录A式（A16）至式（A18），求得各变量仿射方向的增量，并按下式求解迭代步长：

式中：i=1,2，…，r,r为广义不等式约束的维数；ui和wi分别为u和w中的元素；αp为原步长；αd为对偶步长。

步骤6：按式（22）和式（23）计算仿射方向的互补间隙Cgap,af和扰动因子μaf。

式中：Δuaf和Δwaf分别为u和w在仿射方向的增量。

步骤7：按式（24）更新常数项Lu，μ，重新计算Lx′，再次求解附录A式（A16）至式（A18），得到总的牛顿方向，计算相对应的原、对偶步长，最后按照式（25）和式（26）更新所有变量和拉格朗日乘子。

式中：ΔUaf和ΔWaf分别为U和W在仿射方向的增量。

步骤8：令K=K+1，若K≤Kmax，转步骤2继续；否则输出“不收敛”，结束程序。

4 算例分析

本文在MATLAB 2013b平台上实现了最优潮流的各算法编程，并采用Matpower优化软件包[20]中提供的IEEE 30,118,300节点系统以及Polish2 736,3 120节点系统进行了测试分析。表1给出了各个系统的基本参数。

本文程序中涉及的相关参数和初值设置如下：(1)收敛精度ε取10-6;(2)发电机有功出力取上下限的平均值，各节点的电压相角采用平启动，均设为0;(3)松弛变量u取不等式上下限与初值之差的绝对值，拉格朗日乘子y全取1、w全取-1，测试表明松弛变量和拉格朗日乘子的取值对算法影响较大，具体详见附录B。

为了减小计算随机性对分析结果的影响，本文中各算例的CPU平均运行时间为200次运行时间的平均值。由于简化内点算法相比于传统内点算法更为高效，因此，本文所有程序均采用SPCIPM进行求解，且将本文ACOPF和DCOPF程序与Matpower进行对比，验证了本文程序的正确性[20]，具体详见附录C和附录D。

表2给出了各个测试系统基于交流模型、标准直流模型和拟直流模型下最优潮流的最优费用以及与交流模型相比2种直流模型最优费用相对误差的对比分析。与ACOPF的最优费用比较可知，2种直流模型都能够很好地收敛至最优解附近，验证了采用直流模型来近似交流模型的可行性和有效性。DCOPF的最优费用相对误差都保持在3%以内，然而随着测试系统的不断增大，标准直流模型的相对误差也越来越大，对于目前日益庞大和复杂的电力系统而言，DCOPF难以满足大规模电力系统的精度要求。除了IEEE 118节点系统算例以外，PDCOPF最优费用的相对误差都小于1%，相比于DCOPF，计算精度有了很明显的提高，并且系统规模越大，计算精度越高，因此可以很好地应用于现代大规模电力系统最优潮流的在线计算和分析。

在表2的基础上，表3给出了以交流模型运行结果为基准的情况下，标准直流模型和拟直流模型的最优潮流的迭代次数和CPU平均运行时间的对比分析。从表3中迭代次数的对比来看，ACOPF的迭代次数都保持在20次以内，而DCOPF和PDCOPF能在10次内收敛，这是由于交流模型中采用严格的非线性约束数学模型，从而迭代次数明显大于2种直流模型。对比2种直流模型，可以看出PDCOPF的迭代次数比DCOPF减少了约20%。此外，由于采用了具有多项式时间特性的内点算法，3种基于不同模型的最优潮流的迭代次数随测试系统的变化都不是很大。

从表3中的运行时间和时间减少百分比的对比来看，相比于交流模型，DCOPF的运行时间减少了约70%～80%,PDCOPF减少了约60%～70%，可以看出，直流模型比交流模型有着明显的时间优势，尤其是应用于大规模电力系统中。这是由于ACOPF存在非线性约束，运行时间受系统规模的影响较大，而直流模型中由于只包含线性约束，不需要计算各约束的海森矩阵，程序得到很大简化，可以大大减少迭代次数，且运行时间受系统规模的影响较小。虽然PDCOPF的运行时间相对于DCOPF有所增加，但随着系统规模的变大，其相对于ACOPF的时间减少百分比的增加更为明显，因此，在大规模电力系统的应用的优势更为明显。

为了更直观地分析基于不同模型的最优潮流的收敛性能，图1给出了Polish 2 736,3 120节点系统在3种不同模型下最优潮流的迭代收敛过程的比较。可知，内点算法显示出阶段性收敛特性，在收敛至最优解邻域附近时，会进入快速收敛阶段，直至算法完全收敛。基于内点算法在可行域内具有较强的寻优能力，所有模型几乎都能在6～8次迭代之内收敛至最优解邻域，之后进入快速收敛阶段，因此不同模型的优劣决定着内点算法在快速收敛阶段的收敛能力。在快速收敛阶段，基于交流模型的算法需要经过8～10次迭代才能收敛，基于标准直流模型的算法需经过4～5次迭代才能收敛，而基于拟直流模型的算法几乎都在2～3次迭代之内收敛，体现了拟直流模型相对于其他2种模型有着更好的收敛性。

综上，针对现代大规模电力系统而言，ACOPF可以求得精确的数值最优解，但是计算效率偏低，实时计算能力较差；DCOPF比ACOPF有着很好的收敛性和很高的效率，但是计算误差偏大；PDCOPF在满足电力系统在线计算的时间要求的基础上，比DCOPF有着更高的精度，更加符合电力系统在线计算的要求。

5 结论

1）在标准直流模型的基础上，将无功功率作为有功平衡方程的修正量引入有功约束中，建立最优潮流的拟直流模型，大大提高了基于直流模型的最优潮流的计算精度。

2）采用了一种SPCIPM求解拟直流模型的最优潮流，通过优化不等式约束，简化了内点算法的求解过程，便于编程的同时进一步提高了算法的计算效率。

3）通过3个IEEE算例和2个Polish算例的仿真测试，验证了基于SPCIPM的PDCOPF算法的可行性和有效性，其精度和效率可以很好地满足电力系统在线计算和调度的要求。

对偶内点法 篇6

1984年在美国Bell实验室工作的印度学者N.K.Karmarkar[1]提出了“投影尺度法”, 使线性规划出现了自单纯形方法问世以来最重要的一次突破。该算法不仅在理论上被证明是多项式算法, 而且在实际计算中也能与单纯形法相媲美, 尤其在大型问题的应用中, 它显示出能与单纯形法竞争的潜力, 从而打破了单纯形方法一统天下的格局。该方法的出现激起许多学者对内点算法的研究热情。以后许多学者将其算法进行了推广, 用来解决一些非线性规划问题。1989年, R.C.Monteiro和I.Adler[2]给出了用原对偶内点算法求解凸二次规划问题的一个版本, 并证明了其复杂度为O (n3L) , 其中迭代次数为O (nL) 。由于该算法对初始点的要求很严格, 这就给数值实验带来很大的困难。本文我们利用算法的基本思想, 在文献[3]给出线性规划原对偶算法数值实验的基础上, 给出凸二次规划原对偶算法的数值实验, 作为对原文的一个重要补充, 使得内点算法的应用更加广泛。

1 问题提出

凸二次规划原-对偶内点算法的思想是线性规划原-对偶内点算法的一种合理推广[4]。考虑下面的凸二次规划问题 (P) 及其对偶问题 (D) :

这里A是m×n阶的矩阵, 且假设A是行满秩矩阵, Q是n×n阶正定矩阵, b和c分别是m维和n维向量, z是对偶问题中加入的n维松弛变量。

对原对偶问题做出如下假设:

(i) 集合

(ii) 集合, 定义W={ (x, y, z) ;x∈S, (y, z) ∈T}, 和在线性规划中一样, 将对数障碍函数法扩展到凸二次规划中, 通过引入障碍因子μ>0, 我们有下面的问题:

s.t Ax=b, x>0, 其中μ>0, μ称为壁垒参数.则可得到该问题的K-K-T条件如下:

其中, X=diag (x1, x2, …, xn) , Z=diag (z1, z2, …, zn) , e= (1, 1, …, 1) T, 记上述系统的解 (x (μ) , y (μ) , z (μ) ) =ΔW (μ) ,

此时, 原问题与对偶问题的对偶间隙为:

显然当μ※0时, p (x) -d (x) ※0, 因此我们得到下面的定理:

定理1在假设条件 (i) , (ii) 下, 上述系统的解W (μ) , 当μ※0时, x (μ) , (y (μ) , z (μ) ) 分别收敛到原问题 (P) 和对偶问题 (D) 的最优解。

2 算法分析

为了给出算法, 对初始点的选择要求很严格:

原对偶内点算法要求初始点严格满足原可行性和对偶可行性, 而且要求初始点和障碍因子必须满足如下条件:

其中, 是欧式范数, θ是一个很小的正数, 如0.1

算法按下列要求更新障碍因子:

μk+1=μk (1-δ/n) , 其中δ是一个很小的正数, 如0.1。

算法的具体步骤如下:

Step1:选取初始点W 0= (x0, y0, z0) ∈W, 且满足 (*) 式, 给定θ, δ是小的正数, 如0.2, 0.1;给定终止误差ε>0, 置k=0;

Step2:计算对偶间隙 (x (k) ) Tz (k) ≤ε, 则停止迭代, x (k) 为 (P) 的近似解;否则转第三步;

Step3:令μk+1=μk (1-δ/n) , 计算ΔW (k) = (Δx, Δy, Δz) ;

Step4:W (k+1) =W (k) -ΔW (k) , k=k+1, 转第二步;

注1ΔW (k) = (Δx, Δy, Δz) 由用牛顿法解如下方程组确定:

注2关于该算法的收敛性参考文献[2,3,5]。

3 数值实验

例1考察如下的凸二次规划问题及其对偶问题:

给出一组初始解x (1) = (1/2, 9/2, 11/2) T, y (1) = (-15, 50) T, z (1) = (18, 2, 2) T。

取δ=0.1, θ=0.18, 则满足初始解的要求, 用MATLAB 7.0编程, 取ε=0.001, 运行结果如图1所示, 从运行结果可以看出, x越来越趋近最优解 (0, 5, 5) , 从函数图1上可以看出最优目标解趋近235。

例2考虑如下凸二次规划问题及其对偶问题:

给出一组初始解x (1) = (6, 4, 11, 1) T, y (1) = (-3, -3) T, z (1) = (2, 3, 1, 13) T, 取δ=0.1, θ=0.12, 则满足初始解的要求, 用MATLAB 7.0编程, 取ε=0.01, 运行结果如图所示, 过程同例1, 从运行图像上我们可以看出x越来越趋近接近最优解, 从函数图2上可以看出最优目标解趋近。

4 结束语

本文我们初步给出了文献[2]中关于凸二次规划原对偶内点算法的数值例子, 实验结果表明该算法从线性规划推广到凸二次规划后, 仍然具有很好的收敛性和稳定性。

参考文献

[1]Karmarkar N.A new polynomial time algorithm or linear program-ming.Combinatorica, 1984;4 (4) :373—395

[2]Monteiro R D C, Adleri.Interior path following primal-dual algorithm II:convex quadratic programming.Math Programming, 1989;44 (1) :43—66

[3]雍龙泉, 线性规划的原-对偶内点算法数值实验初步.科学技术与工程, 2007;7 (18) :4576—4579

[4]方述成, 普森普拉S.线性优化及扩展, 理论与算法.北京:科学出版社, 1994

对偶内点法 篇7

智能电网的目标是实现电网运行的可靠、安全、经济、高效、环境友好和使用安全。电力系统是一个动态系统,其安全问题的本质必然是动态的[1]。常规的最优潮流(Optimal Power Flow,OPF)难以保证电网运行的动态安全性,不能满足智能电网的要求。暂态稳定约束最优潮流(Transient Stability Constrained Optimal Power Flow,TSCOPF)是一种兼顾安全性与经济性的暂态稳定预防控制实现方法,得到了研究人员的广泛关注[2,3,4,5,6,7,8,9,10,11,12,13,14]。

TSCOPF是在常规OPF中引入一组等式和不等式约束:等式为系统的转子运动方程以及初值方程;不等式为功角稳定性约束。实际上,TSCOPF的约束条件为一组微分代数方程。对微分方程的处理目前主要有两种方法:第一种,将微分方程离散化为一系列的代数方程,直接加到OPF中[2,3,4,5,6,7];第二种,采用约束转换技术将TSCOPF问题化为同规模的OPF问题[8,9,10]。第二种方法虽然可以避免第一种方法中出现规模庞大的问题,但得不到原问题的最优解。现代内点法解优化问题主要的计算量是解修正方程。许多大型的电力系统优化问题都利用稀疏技术及节点优化编号来减少注入元以提高求解速度[15],提升的空间有限。文献[6]在第一种方法的基础上,采用预测校正内点法与简约空间技术相结合的算法对TSCOPF问题进行求解,缩减了求解修正方程所耗用的时间,有效地缓解了上述问题。

简约空间技术是一种高效的线性方程求解器,最初被用于化学工程及流程优化等领域[12,13]。简约空间技术按照一定的分解策略对修正方程进行分解,通过缩小修正方程的规模,减少解修正方程所耗用的时间,从而达到提高计算效率、降低内存损耗的目的。简约空间技术适合解自由度小的问题,其计算复杂度主要取决于问题自由度的大小,即控制变量与总原始变量的比值的大小。离散化后的TSCOPF问题虽然规模庞大,但其控制变量的个数相对很小,且不随故障重数、仿真时长和仿真步长的改变而改变,因此适于用简约空间技术求解。

转子运动方程包含了电磁功率方程,大部分文献中的电磁功率方程采用混合极坐标[2,3,4,5,6,7,9,11,12]。与直角坐标及混合极坐标形式相比,极坐标形式的模型最为简洁,便于记忆和编程。本文在文献[6]的基础上,建立了新的基于极坐标形式的TSCOPF模型,即潮流方程及电磁功率方程均采用极坐标。统一的坐标系,降低了程序的编制复杂度。采用了简约空间内点法对该模型进行求解。测试结果表明,该模型简洁可靠,所提算法计算效率高,占用内存少,适于求解大规模的TSCOPF问题。

1 极坐标形式的TSCOPF模型

OPF模型在数学上可归纳为如式(1)所示。

其中:x∈Rn为变量;f(x):Rn→R为目标函数;h(x):Rn→Rm为等式约束;g(x):Rn→Rc为不等式约束;分别为不等式约束的上下限。TSCOPF模型的差异主要体现在电磁功率方程及潮流方程的坐标形式的选取上。

1.1 三种坐标形式比较

以有功功率方程为例,目前混合极坐标、极坐标及直角坐标形式的潮流方程都已被采用,如式(2)~式(4)所示,但电磁功率方程还只是局限于混合极坐标,如式(5)所示。

其中:i∈SN,SN为所有节点集合;SG为发电机节点的集合;PGi、PDi为节点i处有功源输出功率、有功负荷;Ui、θi为节点i电压幅值、相角;Gij、Bij为系统导纳矩阵的实部、虚部;Yij、αij为节点导纳矩阵的幅值、相角;ei、fi为节点i电压值的实部、虚部;Ptei、δit为发电机i在t时刻的电磁功率、转子角度;Ei为发电机i的暂态电抗后电势;Geij、Beij为仅保留发电机节点的既约导纳阵实部、虚部。

文献[6]采用式(4)与式(5)组合模型一,文献[3]采用式(2)与式(5)组合模型二,文献[4]采用式(3)与式(5)组合模型三。组合模型一、模型三中出现两种坐标系,使得模型复杂,后续的公式推导及编程的工作量成倍数增加,且出错率增大。比较式(2)、式(3)与式(4)可以发现,混合极坐标形式的复杂度约为极坐标形式的两倍,其计算量远大于极坐标形式。虽然直角坐标模型的海森矩阵为常数阵,但不如极坐标模型简洁。基于此,本文提出基于完全极坐标形式的TSCOPF模型,即潮流方程及电磁功率方程均采用极坐标。极坐标形式的电磁功率方程如式(6)所示。

其中:Yeij、αeij为仅保留发电机节点的既约导纳矩阵的幅值、相角;其他变量定义与前相同。

1.2 极坐标形式的数学模型

本文考虑经典模型:发电机用暂态电抗Xd′和暂态电抗后电势E′相串联表示,且设E′恒定;发电机输入的机械功率Pm恒定;负荷用恒定阻抗表示。取系统网损最小为目标函数。极坐标形式的TSCOPF模型描述如下:

其中:,ST为积分时段的集合;SR为无功源点的集合;QRi、QDi为节点i处无功源输出功率、无功负荷;为发电机i在t时刻的角速度;Δt为积分步长;iT为发电机i的转动惯量;Xd′i为发电机i的暂态电抗;分别为有功源出力上下限、无功源出力上下限和节点电压幅值上下限;为发电机转子相对摇摆角度上下限,取为±100°;δ0COI为稳态时惯性中心角度;δtCOI为预想故障下各时段的惯性中心角度;ωN为角速度基准值;其他变量定义与前相同。

比较各模型可见,本文模型最为简洁,便于记忆及编程。统一的坐标系,使后续的梯度矩阵、海森矩阵的公式推导及编程实现的工作量成倍数减少,有效地降低了出错率,提高了计算效率。

2 简约空间内点法求解步骤

本节首先介绍简约空间技术的推导过程并分析其计算复杂性,然后阐述简约空间内点法解TSCOPF问题的具体步骤。

2.1 简约空间技术公式推导

线性系统问题的标准格式为

其中:D∈Rn×n;A∈Rm×n;Lx∈Rn;Ly∈Rm。

将原变量x分解为两个子空间DY和DZ,令

由式(17)即可得Δy。

2.2 计算复杂性分析

内点法既约修正方程即为如式(7)所示的线性系统。用符号W(·)表示式(7)中系数矩阵,则W(·)∈R(n+m)×(n+m)。W(·)为高度稀疏对称阵,常规的处理方法为采用LDLT分解直接求解。该部分的计算量消耗了整个程序的绝大部分计算时间,其阶数直接影响到内点法的计算效率。

简约空间技术的主要计算量集中于求解式(10)、式(14)、式(15)和式(17),其计算复杂度主要取决于问题自由度的大小。当系统自由度很小,即n-m/n<<1时,系数矩阵C∈Rm×m,B∈R(n-m)×(n-m),修正方程的阶数被大大缩小,能有效地减少解修正方程所耗用的时间。且观察式(10)、式(15)、式(17)的系数矩阵发现,对该三式只须进行一次系数矩阵的分解,可采用前代、回代对其他两式进行求解。

2.3 简约空间内点法求解步骤

原始-对偶内点算法[18]在电力系统中得到了很好的应用。文献[3]亦证明了原始-对偶内点法求解TSCOPF问题的有效性。基于简约空间内点法的TSCOPF问题的求解步骤如下。

步骤1:读入系统数据及故障信息,置kmax=50,中心参数取σ=0.1,计算精度为ε=10-5。

步骤2:常规潮流计算,得到电压初值。

步骤3:根据常规潮流结果计算各负荷等值阻抗,结合预想故障信息,消去网络节点,形成故障中及故障切除后仅含发电机节点的既约导纳阵。

步骤4:设定k=0。选择初值l>0,u>0,z>0,w<0,y=0。

步骤5:判断k≤kmax是否成立,如果不成立,转步骤13。

步骤6:计算互补间隙,如果CGap<ε,则输出最优结果,停机。

步骤7:计算扰动因子

步骤8:计算KKT条件中的Ly、Lw、Lz、Ll、Lu,目标函数、等式及不等式约束的雅克比矩阵和海森矩阵∇f(x)、∇h(x)、∇g(x)、∇2f(x)、∇2h(x)、确良∇2g(x),构造既约修正方程中的H′、Lx。既约修正方程为

其中:

步骤9:简约空间技术解修正方程。

步骤10:计算

步骤11:求最大原始和对偶步长。

步骤12:更新变量

k=k+1,转步骤5。

步骤13:提示“计算不收敛!”,停机。

3 测试结果与讨论

本节以表1所示的五个测试系统为例,验证上述简约空间内点法解极坐标形式的TSCOPF问题的有效性。从表1可以看出,各系统TSCOPF问题的自由度都很小。各系统预想故障信息如表2所示。目前只考虑单重故障。故障节点的故障类型为三相短路。采用32位Matlab R2010b。计算环境为Dual-Core AMD 2.2 GHz,2 GB内存。

互补间隙是衡量系统最优性的重要指标。各系统的互补间隙随迭代次数的变化曲线如图1所示。从图1可以看出,各测试系统的对偶间隙单调收敛至零,该算法具有良好的收敛性。

以CEPRI-22、NE-39和IEEE-118系统为例,发电机转子相对摇摆曲线如图2所示。所有发电机的转子相对角度均在-100°~100°范围内,满足了在仿真时间内各发电机不失稳的要求。

常规原始-对偶内点法与简约空间内点法的测试结果比较如表3所示。从表3可以看出,简约空间内点法的CPU总计算时间及解修正方程所耗用的时间均远小于常规原始-对偶内点法。

若定义简约空间内点法的CPU计算时间与常规原始-对偶内点法的CPU计算时间的比值为加速比,则测试系统的加速比如图3所示。可见,加速比曲线随着系统规模的增大呈单调递增趋势,系统规模越大,简约空间内点法加速越明显。

图4给出了常规原始-对偶内点法和简约空间内点法解修正方程的时间占CPU总计算时间的百分比。两条百分比曲线均随着系统规模的增大呈单调递增趋势,可见,系统规模越大,内点法中修正方程的求解问题越严重。由简约空间内点法的百分比曲线低于常规原始-对偶内点法的百分比曲线可以看出,简约空间内点法有效地缩减了解修正方程所耗用的时间,提高了计算效率。

4 结论

本文提出了极坐标形式的TSCOPF模型,该模型将潮流方程及转子运动方程统一到极坐标形式下,比其他形式的模型更简洁明了,便于记忆及编程。采用了简约空间内点法进行求解,这一方法通过缩小修正方程的规模,减少解修正方程所耗用的时间,达到了提高计算效率、降低内存损耗的目的。