基于凸包理论的含风电电力系统负荷恢复方案优化

2024-03-29伊昆明

电力系统自动化 2024年5期

伊昆明，孙磊，丁江，丁明

（1.合肥工业大学电气与自动化工程学院,安徽省合肥市 230009；2.新能源利用与节能安徽省重点实验室（合肥工业大学）,安徽省合肥市 230009；3.国网南昌供电公司,江西省南昌市 330000）

0 引言

近年来,自然灾害和网络攻击等极端事件导致电力系统大规模停电事故频发,给国民经济带来严重的负面影响［1］。合理有效的负荷恢复方案有助于加快电力系统整体恢复进程,大幅降低因大停电事故导致的负面影响［2］。以风电为代表的可再生能源接入电力系统将有利于加快系统恢复速度［3］,但其出力的不确定性使得系统在恢复过程中存在大停电风险［4］。因此,研究含风电的电力系统负荷恢复策略具有重要意义。

现有负荷恢复模型通常考虑频率、电压、发电机出力和潮流方程等约束条件,以最大化加权负荷恢复量或最小化停电损失为目标函数［5-8］。文献［9］以最大化负荷恢复量为目标,建立了基于交流（AC）潮流的负荷恢复模型,但无法保证模型在多项式时间内精确求解。为解决考虑AC 潮流的负荷恢复模型求解效率低的问题,已有专家学者展开了对AC 潮流方程的线性近似和凸松弛的研究。线性近似方法主要对方程中的非线性项做近似处理,将非凸潮流模型转化为线性的潮流模型［10］,如直流潮流模型。文献［11-12］基于直流潮流提出了考虑动态频率安全约束的负荷恢复模型。然而,直流潮流模型忽略了无功功率和网损,导致所得解误差较大。文献［13-14］将正余弦函数和电压的平方项描述为关于相角和电压的线性函数,建立了基于混合整数线性规划的滚动负荷恢复模型。但线性化的模型所得到的解未必是原模型的可行解［15］。凸松弛方法主要包括半正定松弛（semidefinite relaxation,SDR）、二阶锥松弛（second-order cone relaxation,SOCR）和二次凸松弛（quadratic convex relaxation,QCR）,其中,SOCR 因其求解速度较快而被广泛应用在电力系统诸多领域［16-17］。文献［18］提出了计及通信系统失效的信息物理协同恢复策略,并构建了基于SOCR 的恢复模型,提高了模型的求解效率。文献［19］建立了以最大化弹性指标为目标函数的故障恢复模型,并将所提出的模型描述为混合整数二阶锥优化模型。SOCR 具有快速求得全局最优解的优势,但由于松弛误差较大,其得到的最优解未必在原模型可行域内［15］。采用上述文献的方法处理负荷恢复模型中的AC 潮流方程,虽然在一定程度上提高了模型的求解效率,但无法保证解的精度。

凸包理论在处理强非线性电力系统AC 潮流方程方面得到了一定应用。文献［20］提出了单变量平方、双变量乘积等非线性项的凸包方程。在此基础上,文献［21-22］将非线性项的凸包方程应用到最优潮流中,算例验证了凸包方法能够减小松弛解与标准模型解的最优性差距。文献［23］提出了将潮流方程的极坐标形式转化为复功率指数形式的坐标变换方法,所提出的方法有助于提高凸包紧密性。文献［24］提出了基于正交和组合的方法,有效减少了AC 潮流SDR 的误差。极坐标下的AC 潮流方程具有高度非线性,由递归方法得到的三变量乘积项的凸表达式无法保证是三变量乘积项的凸包［25］,其带来的松弛误差较大。凸包方法的效果在一定程度上取决于凸包参数,但目前尚且没有精确的凸包参数的确定方法,通过调整凸包参数的取值范围,可缩小可行域,进而减少因引入凸包带来的误差。

针对上述问题,本文首先针对AC 潮流方程的非凸性,建立了基于直角坐标的负荷恢复二阶锥优化模型。其次,采用条件风险价值（conditional value-at-risk,CVaR）理论处理风电的不确定性。然后,基于凸包理论对引入的松弛变量进行凸化,构建基于凸包理论的混合整数二阶锥优化模型,优化负荷恢复方案；提出凸包参数的动态迭代方法,缩小可行域范围,减少松弛误差。最后,采用IEEE 39节点系统和改进的广东电力系统,验证所提模型的有效性。

1 凸包理论

凸包本质上是一个凸多边形或多面体,通过多个凸函数和凹函数包围一个非线性方程。其中,凸函数是不大于非线性方程全局解的最高凸函数,凹函数是不小于非线性方程全局解的最低凹函数,由这些凸函数和凹函数构成的集合称为凸包［20］。凸包的意义在于能够对非线性约束进行凸松弛,对约束中非线性项凸化。电力系统的AC 潮流方程包含三角函数、变量相乘等非线性项,导致模型求解困难且无法保证全局最优。将SOCR、QCR 等凸松弛技术引入电力潮流方程中已成为研究热点,因其能够加快求解速度且在一定松弛条件下能够求得全局最优解而被广泛应用。常规SOCR 方法虽能提高模型求解效率,但当其松弛误差较大时,所得结果并不精确,且仅能得到原问题最优值的一个上界或下界［15］。基于凸包理论,通过动态改变包络边界缩小可行域,提高松弛紧密性,减少松弛后的最优解与原模型最优解的误差,从而提高模型求解的精度,且运用凸包理论得到的解与原问题的解相近甚至相同。

确定函数在其定义域内的凸多面体（convex polyhedron）有利于构造该函数的凸包方程,这是因为凸多面体包络（convex polyhedron envelope）总可由有限个仿射函数表示,且该包络为所有仿射函数在其定义域内的最大值［26］。仿射函数具有F(x)=ATx+b形式,其中,A和x为n维向量,b为参数。具有凸多面体的函数,其凸包一定经过多面体上的顶点,本文通过可行域的顶点获取超平面进而得到凸包方程。虽然文献［23］提出了定义域内分区域构造凸包的方法,即将定义域分成若干子定义域,在子定义域内由顶点获取超平面,但会增加模型中布尔变量和约束的个数,导致模型求解难度增加。因此,确定一个具有凸多面体的函数的凸包方程,需要找到该函数在定义域内的顶点,例如,节点i电压Vi的平方 (∀i∈ΩB) 在定义域内的顶点为(Vmin,VminVmin)和(Vmax,VmaxVmax),其凸包表达式见式（1）和式（2）；双变量乘积项ViVj(∀i,j∈ΩB)在定义域内的顶点为(Vmin,Vmin,VminVmin)、（Vmax,Vmax,VmaxVmax）、（Vmin,Vmax,VminVmax）和（Vmax,Vmin,VmaxVmin）,则该双变量乘积项ViVj的凸包表达式如式（3）—式（6）所示［20-21］。

式中：ΩB为节点集合；Vmin和Vmax分别为节点电压的最小和最大值；uii和uij为表示凸包而引入的辅助变量。

由式（1）和式（2）构成的松弛域见附录A 图A1,由式（3）—式（6）构成的松弛域见附录A 图A2。附录A 图A1 展示了V2i的凸包松弛域,如阴影区域所示；附录A 图A2 展示了ViVj的凸包松弛域的一个截面,其中,阴影区域表示的凸松弛域由4 条直线确定其边界。

文献［26-27］证明了函数ψ(y1,y2)=y1y2在定义域R2中的所有多面体上均有一个凸多面体包络,且函数在该凸多面体是边凹（edge concave）［26］的。因此,其凸包方程可以直接由凸多面体顶点确定［28］。设（k′为边的条数）为多面体P所有边的集合,SP为多面体P所有点的集合,函数G在多面体P上二阶连续可微,HG(p)为函数G在点p处的Hessian 矩阵,若对于D中任意边向量dli都有,其中,p∈SP,则称函数G在该多面体上是边凹的［26,29］。如果函数G在它的一个多面体上均是边凹的,则该多面体为凸多面体且有一个凸多面体包络。如果一个函数φ(m1,m2,…,mn)对于任一变量mi是单独线性的（其他变量设定为参数,n为变量的维数）,即mi为变量时,mj(j≠i)均为参数,则函数φ(m1,m2,…,mn)称为多线性函数（multilinear function）。已有文献证明一个多线性函数在其定义域内有凸多面体包络,且函数在多面体上是边凹的［26,29］。

由上述理论可知,函数Θ(y1,y2,y3,y4)=y1y2+y3y4是由形如函数ψ(y1,y2)=y1y2线性叠加而成的,且Θ(y1,y2,y3,y4)对于每一个变量都是单独线性的,故其为多线性函数。Θ(y1,y2,y3,y4)在定义域上的凸包方程可以直接由凸多面体顶点确定。如果函数ψ和函数ϕ的凸包均是凸多面体,则函数ψ的凸包与函数ϕ的凸包的Minkowski 和,即conv(ψ)+conv(ϕ)也是凸多面体［30］。由文献［26］可知两个函数的凸包之和应包含在这两个函数和的凸包内,如式（7）所示。

式中：ψ和ϕ均为表示含双变量乘积项的函数；conv(⋅)为函数的凸包。当且仅当函数ψ的凸包与函数ϕ的凸包的Minkowski 和以及函数Θ的凸包均为凸多面体时,式（7）中的等号成立。

式（1）和式（2）可构建平方项函数的凸包；式（3）—式（6）可构建双变量乘积项的凸包,结合式（7）及其等号成立条件,可构建函数和的凸包,将在2.2节中具体介绍。附录A 图A3 展示了函数的凸包松弛域,为附录A 图A3 中蓝色平面与红色曲面之间的区域。其中,ei,t和fi,t分别为第t时步电力系统节点i的电压实部和电压虚部。

2 基于凸包理论的负荷恢复模型

在实际恢复过程中,负荷的投入与否由断路器的开合状态决定,而一条母线通常有多条馈线。因此,节点负荷的恢复量不可能连续变化。根据电力系统实际运行工况,母线负荷由若干离散的负荷点组成。构建的负荷恢复模型中决策变量分别为每个时步负荷点的恢复状态和发电机出力等。

2.1 目标函数

负荷恢复阶段的主要目标是尽快恢复电力负荷,本文所提出的负荷恢复模型的目标函数之一为最小化停电损失。风电作为一种发展较为成熟的清洁能源,已广泛接入电力系统,但风电机组出力的不确定性影响系统稳定运行和恢复能力。本文引入节点可容许风电功率范围的概念［31-32］,并结合CVaR理论计算因风电出力的不确定性带来的风险损失。因此,本文模型的另一个目标函数为最小化风电功率缺额和盈余风险损失值。综上所述,本文提出的目标函数可由式（8）—式（10）表示。

本文模型中的变量和参数如附录B 表B1 和表B2 所示。去除式（8）中的常数项,目标函数可表示为：

式（9）和式（10）包含非线性积分项,对其线性化,分别如式（12）、式（13）和式（14）、式（15）所示。

2.2 约束条件

1）节点平衡方程和潮流方程。

采用文献［33］提出的二阶锥优化的方法处理潮流方程,引入辅助变量wii,t、cij,t、sij,t,将AC 潮流方程描述为二阶锥形式。引入的辅助变量分别如式（16）—式（18）所示。

式中：ΩL为输电线路集合。

节点平衡方程和潮流方程可描述为：

式中：ΩGi为节点i的火电机组集合；和分别为节点i的并联电导和并联电纳；和分别为节点导纳矩阵中第i行第j列元素对应的电导和电纳；为线路ij的充电电纳；Q为节点i的负荷d的无功功率；P和Q分别为节点i的火电机组g在第t时步的有功功率和无功功率；Pij,t和Qij,t分别为第t时步线路ij的有功和无功功率。式（19）和式（20）为节点平衡方程。式（21）和式（22）为输电线路潮流方程。

2）节点电压约束见式（23）,可约束电压范围。

3）机组出力与爬坡约束见式（24）—式（26）。

4）电力线路传输功率约束见式（27）,表示线路上流过的功率应在给定的阈值范围内。

式（27）可线性化为：

5）平衡节点电压约束见式（32）,表示平衡节点的电压虚部为0。

式中：∀i：ref 表示可以任意指定节点i为参考节点。

6）最大负荷恢复量约束见式（33）和式（34）。

式中：ΔPmaxt为第t时步允许的最大负荷恢复量；ρ为最大允许负荷恢复量与总的发电机有功出力的比值,其取值通常为5%［9］,对其取值的讨论可参考文献［34］。式（33）表示在第t时步负荷恢复量不能超过允许的最大恢复量。式（34）表示在第t时步最大负荷恢复量受到t-1 时步总的发电机有功出力限制。

7）旋转备用约束见式（35）—式（40）。

式中：μ为P的最大值与的比值；R和R分别为火电机组g在第t时步的向上和向下备用容量；ΔPDt为第t时步系统的负荷恢复量。式（35）表示机组向上备用容量提供给负荷恢复和风电功率缺额。式（36）表示机组向下备用容量提供给风电功率盈余。式（37）和式（38）表示功率可容许范围下限和上限均应在阈值范围内。式（39）和式（40）表示机组当前时步所能提供的向上和向下备用容量。式（39）和（40）可采用大M 法线性化,此处不再赘述。

8）负荷恢复状态约束见式（41）,表示负荷点一旦恢复,其恢复状态将保持不变。

9）辅助变量的凸包约束见式（42）—式（44）。

根据辅助变量wii,t、cij,t和sij,t的定义,可知三者满足：

式（42）为松弛的二阶锥表达式。式（43）和式（44）分别为辅助变量cij,t和sij,t应满足的约束关系。

为减少基于SOCR 的负荷恢复模型带来的误差,基于第1 章介绍的凸包理论,对辅助变量wii,t、cij,t和sij,t进行凸化,如式（45）—式（62）所示。

综上所述,基于凸包理论的负荷恢复模型可描述为：

3 求解步骤

为了提高负荷恢复模型的求解精度,本文提出直角坐标系下基于凸包理论的负荷恢复二阶锥优化模型,并提出约束收紧的迭代方法动态修改模型中的参数。基于参数迭代的电力系统负荷恢复模型的求解流程如图1 所示,具体步骤描述如下：

图1 基于参数迭代的电力系统负荷恢复模型的求解流程图Fig.1 Flow chart of solution for load restoration model of power system based on parameter iteration

步骤1：输入参数松弛容差δ、松弛因子ε(k)、调节系数γ和最大迭代次数K,初始化凸包参数、,令k=1。

步骤2：求解式（63）所示模型。

步骤3：根据求解器输出的求解状态,判断是否得到最优解,若是,执行步骤4,否则转至步骤7。

步骤4：根据式（64）,计算每条线路的相对松弛误差ηij,t,判断所有线路的相对松弛误差是否均不大于允许误差或迭代次数是否大于K+1,若是,则求解流程结束并输出负荷恢复方案,否则转至步骤5。

步骤5：将步骤2 中得到的ei,t和fi,t分别赋值给中间变量e和f。

步骤6：基于松弛因子ε(k)更新凸包参数；令k=k+1,转至步骤2。

步骤7：更新第k次迭代的ε(k),转至步骤6。

4 算例分析

本章以IEEE 10 机39 节点系统和改进的广东电力系统为例来说明所提出方法的有效性。本文建立的模型为二阶锥优化模型,程序基于AMPL 软件实现,通过调用Gurobi 9.5.2 求解器对所提出的模型进行求解。运行环境是处理器为Intel Core i7-10700、内存为16 GB 的台式电脑。

4.1 IEEE 10 机39 节点系统

IEEE 10 机39 节点系统的负荷参数和机组参数分别见附录C 表C1 和表C2,机组启动功率可详见文献［35］,其中,待恢复负荷总量为2.77 GW。离散化的恢复时间间隔为10 min,总恢复时长设置为140 min,ρ取值为5%。节点3、5、14 和16 处有4 个风电场,各风电场的预测出力曲线见附录C 图C1。c0、c1、c2分别设置为1 4000、500、50［30］,μ取值为1。

4.1.1 负荷恢复结果

附录C 表C3 所示为求解本文模型得到的最优负荷恢复方案。由表C3 可以看出,在得到的恢复方案中,恢复所有负荷所需时间为130 min,且各时步负荷恢复量呈现递增趋势。然而,相比第110 min时的负荷恢复量,第120 min 和130 min 时的负荷恢复量有所减少,这是因为负荷恢复后期,待恢复负荷较少,且由于第120 min 和130 min 时风电缺额风险减小会导致负荷恢复量有所减小。

不同风电场在各时步的功率可容许范围的上下限值如附录C 图C2 所示。由图C2 可以看出,所有风电场在各时步的功率可容许范围上下限的差值均大于90 MW,根据不同风电场在各时步的功率预测值和预测误差得到其概率密度函数,并结合式（65）可计算得到风电实际出力在其功率可容许范围的概率。在本算例中,针对所有时步风电实际出力在功率可容许范围的概率均不小于81%。该概率值相对较小,这是因为本算例中设置μ为1,通过减少μ的值,可增大可容许范围,相应的概率也会增加。

4.1.2 风电出力对负荷恢复影响

为了分析风电出力对负荷恢复结果的影响,根据风电是否参与、是否考虑风险以及μ的不同取值设置5 种不同场景,见附录C 表C4。场景1：允许风电参与,但不考虑风险,风电出力为预测值；场景2、3、4：均允许风电参与并考虑风险,区别在于μ的取值,分别设置为1.00、0.85、0.70；场景5：无风电参与。

不同场景的各时步总负荷恢复量对比如图2 所示。由图2 可以看出：第130 min 时,场景1 和场景2的总负荷恢复量达到2 770 MW,负荷恢复过程结束；在第140 min 时,场景3、4、5 完成负荷恢复。相比于场景5,场景1 完成负荷恢复所需的时间较短,进而说明风电的参与能够加快负荷恢复进程。相比于场景2,场景3、4 的μ值有所减少,导致负荷恢复进程延缓。场景2、3、4 的各时步负荷恢复量和容许风电缺额量对比如附录C 图C3 所示。容许风电缺额量是指风电预测值与可容许范围下限之差,其值越大,则由风电缺额带来的风险越小。由附录C 图C3 可以看出：第80 min 至120 min 时段,针对负荷恢复量大小满足场景2>场景3>场景4；针对容许风电缺额量大小满足场景4>场景3>场景2。结合图2 和附录C 图C3 可知：减小μ值能够增大容许风电缺额量,进而降低风电不确定性带来的风险和单个时步负荷恢复量,导致总负荷恢复量有所减少,延缓了负荷恢复进程。

图2 不同场景的各时步总负荷恢复量Fig.2 Total restored load capacity at each time step in different scenarios

4.1.3 本文模型的有效性验证

为了说明本文提出直角坐标系下基于凸包的负荷恢复模型的有效性,将本文所提出的模型与文献［9］和文献［36］的模型进行对比。文献［9］的模型为基于标准AC 潮流方程的负荷恢复模型；文献［36］的模型为基于二阶锥优化的负荷恢复模型。

表1 为求解不同模型得到的目标函数值和求解时间。

表1 求解不同模型得到的目标函数值和计算时间Table 1 Objective function value and computational time by solving different models

由表1 可以看出,求解文献［36］的模型得到的目标函数值最大,求解考虑参数迭代的本文模型得到的目标函数值次之,求解文献［9］的模型得到的目标函数值最小。由此可知,相比于文献［36］的模型,采用本文提出的模型得到的解与文献［9］的模型得到的解的误差更小,其中,误差由式（66）计算。与考虑参数迭代的场景相比,求解不考虑参数迭代的本文模型得到的解的误差更大。由此可知,采用基于参数迭代的方法有助于减少松弛误差。需要指出的是,本文模型虽然牺牲了一定的求解时间,但提高了模型的求解精度,且本文模型的求解时间在可接受范围内。

式中：ε′为相对误差；SAC为求解文献［9］的模型得到的最优解；Srelax为求解松弛模型的最优解,在本算例中,松弛模型包括文献［36］的模型和本文模型。

求解不同模型得到的最后一个恢复时步的线路相对松弛误差如图3 所示。由图3 可以看出,求解文献［36］的模型得到的恢复方案中,60%线路的相对松弛误差大于0.01%；而求解本文模型得到的恢复方案中,所有线路的相对松弛误差均小于0.005%。因此,本文提出的模型能够减小因SOCR 带来的线路相对松弛误差,提高松弛紧密度。

图3 求解不同模型得到的最后一个恢复时步的输电线路相对松弛误差Fig.3 Relative relaxation errors of transmission lines at the last time step by solving different models

求解不同模型得到的最后一个恢复时步的部分线路的有功功率如附录C 表C5 所示。由表C5 可以看出,求解文献［36］模型得到的线路有功功率最大相对误差为10.9%,而求解本文模型得到的线路有功功率的最大相对误差仅为1.1%,满足工程要求。因此,相比于文献［36］的模型,本文提出的模型能够得到更接近求解文献［9］的模型得到的全局最优解。

4.2 改进的广东电力系统

为了进一步说明本文所提出模型的有效性,本节采用改进的广东电力系统对其进行验证。改进后的广东实际电力系统包含339 个节点、611 条线路、51 台火电机组。其中,待恢复负荷为9.68 GW；10 个风电场分别位于节点5、14、45、49、124、126、127、149、157 和158,其预测功率见附录C 图C4。每个恢复时步步长为15 min,总恢复时长设置为210 min,ρ取值为2%。

不同场景的各时步总负荷恢复量如附录C 图C5 所示。由图C5 可以看出,第180 min 时,场景1和场景2 的总负荷恢复量为9.68 GW,负荷恢复过程结束；场景3、4、5 在第195 min 完成负荷恢复。相比于场景5,场景1 完成负荷恢复所需的时间较短,说明风电的参与能够加快负荷恢复进程。相比于场景2,场景3 和场景4 的μ值有所减少,导致负荷恢复进程延缓。

求解不同模型得到的目标函数和求解时间见附录C 表C6。由表C6 可以看出,求解文献［36］模型得到的目标函数值为588 006万美元,求解考虑参数迭代的本文模型得到的目标函数值为587 937 万美元,求解文献［9］模型得到的目标函数值为587 901万美元。由此可知,相比于文献［36］模型,采用本文所提模型得到的解更接近于基于AC 潮流方程的解。求解不考虑参数迭代的本文模型得到的目标函数值为588 001 万美元,其误差比基于参数迭代模型得到的解的误差大。由此可知,采用基于参数迭代的方法有助于减少松弛误差。

不同风电装机占比场景的各时步总负荷恢复量如图4 所示。风电装机占比为风电场装机容量与系统总装机容量的比值,在该系统中除风电外的其他机组装机容量不变,通过改变风电装机容量,得到不同风电装机占比场景,进而分析不同风电装机占比对负荷恢复方案的影响。由图4 可以看出,在风电比例分别为0%、4%、8%、12%、16%和20%时,完成负荷恢复所需的时间分别为195、180、180、165、165、150 min。随着电力系统中风电比例的增加,各时步总负荷恢复量呈递增趋势,完成负荷恢复所需的时间呈递减趋势。由此可知,在电力系统中接入风电将有助于加快负荷恢复进程。