主动配电网二阶段鲁棒日前调度优化模型及其求解方法

时间：TIME\@"yyyy'年'M'月'd'日'"2022年3月29日学海无涯页码：第1-页共1页主动配电网二阶段鲁棒日前调度优化模型及其求解方法随着电力体制改革的不断推进和间歇性分布式电源（distributiongeneration，DG）渗透率的不断增大，配电网需要面对愈加不确定的能源环境。配电网是智能电网实现精益化电能分配和保证电网稳定运行的重要环节，合理制定配电网日前调度计划将是电力系统可持续发展战略的重中之重。

为了应对配电网日前调度计划问题中的不确定因素，目前学者主要提出以下几种方法：随机优化、机会约束算法和鲁棒优化等。随机优化通常采用以文献［1-3］为代表的概率密度分布模型和概率场景模型来描述分布式电源和负荷的不确定性，再通过蒙特卡洛或场景分析等方法将不确定问题转化为确定性问题。以文献［4］为代表的机会约束算法延续了概率场景模型的思路，额外引入置信度指标，调度优化模型中的约束条件只须在一定置信水平成立即可。

然而以上两种方法都需要大量的历史数据拟合出概率分布模型，且随着求解维度的扩增，求解时间也呈指数级增长。相比之下，鲁棒优化方法只需要获取通过少量历史数据确定的区间分布情况，具有更高的工程应用价值［5］，且以文献［6］为代表的可调鲁棒优化算法，使得模型能在鲁棒性和经济性之间做出权衡。

鲁棒优化在配电网优化方面逐渐得到重视，近年来二阶段鲁棒优化算法成为主要研究热点，例如文献［7］提出考虑风电和需求响应的机组组合二阶段鲁棒优化方法，以及文献［8］提出含电动汽车充电站的主动配电网二阶段鲁棒规划模型，但是二阶段鲁棒优化算法在配电网日前计划调度方面的研究刚步入初始阶段。文献［9］通过量化调度者多元化风险偏好，可适当调整优化方案的保守性，进而利用主‐子问题迭代算法求解二阶段鲁棒模型的最优值，但该模型过度简化，未建立不确定变量与主‐子问题决策变量之间的约束关系式，因此忽略了子问题在强对偶过程中形成的双线性非凸因素。文献［10］引入big-M法对子问题经对偶转化后的双线性项进行处理，但同时给子问题引入了整数变量，有悖于列约束生成算法中子问题只能是线性规划的必要条件，将无法保证子问题返给主问题的是最优割平面。文献［11］将个别双线性项中的不确定变量取值为其波动区间的边界值，符合鲁棒优化中“最恶劣”场景的定义，但是该最优潮流部分仅仅考虑了有功潮流，应用于求解配电网日前调度时会存在较大的误差。

针对以上不足以及考虑到二阶段鲁棒日前调度模型为min-max-min结构的混合整数规划问题，本文借鉴文献［12］提出的基于Benders分解法和外层逼近法（outerapproximation，OA）相结合的求解方法，首先通过外层逼近法处理子问题对偶过程中产生的双线性项，再用Benders分解算法对主‐子问题进行迭代求解，同时引入文献［13］中的二阶锥模型和其对偶转化形式。

此外，国内目前很少在配电网二阶段鲁棒调度模型中考虑电价波动，电价由政府调控，用电弹性需求较低，因此学者普遍将主网的分时电价作为已知参数。文献［14］在峰谷分时电价的基础上制定实时电价，引导用户响应电价的变化以调整用电需求。但是随着新一轮电力体制的改革，研究配电网调度优化在电力市场背景下的挑战与应对策略将至关重要，尤其是电价波动的影响。在自由化的电力市场中，电价由供需双方决定，并且由于电力商品具有不可贮存、实时平衡、输送约束等特性使得电价存在波动性［15］。文献［16］虽然提出考虑电价波动性的鲁棒优化模型，但其研究主题是以配电网的网损成本为目标函数的最优潮流算法，属于配网调度的小范畴。因此本文在电力市场背景下提出考虑电价波动的配电网日前调度模型。

1配电网日前调度优化模型

图1为典型的配电网日前调度优化结构图。本文建立的配电网日前调度优化模型中，主网电价信息、负荷预测信息和DG预测信息为调度中心的决策依据，主网购电量、储能充放控制策略、可控分布式电源发电量和无功补偿量为调度中心的决策变量，以配电网的总运行成本最小为模型寻优目标，同时调控策略要确保配电网侧和可控资源侧满足安全运行的约束条件。

图1配电网日前调度优化结构图

Fig.1Day-aheaddispatchingoptimizationstructurediagramofdistributionnetwork

1.1目标函数

模型的目标函数为配电网调度成本函数，该函数由配电网侧产生成本和可控资源侧成本两部分构成。配电网侧部分涵盖主网功率购电费用和网损费用，可控资源侧涵盖可控式储能单元运行费用、传统分布式发电设备的运行燃料费用和新能源出力的削减费用，其中传统分布式发电设备可采用微型燃气轮机，其发电成本可由线性函数表示［11］。调度优化模型的目标函数F具体表示如下：

F=min(CM+CLoss+CBat+CG+CDG)

(1)

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪CM=∑t∈ΩT∑i∈ΩMcMtPMi,tΔtCLoss=∑t∈ΩT∑ij∈ΩlinecLossI2ij,trijΔtCBat=∑t∈ΩT∑i∈ΩCcdisPdisi,tΔt+cchPchi,tΔtCG=∑t∈ΩT∑i∈ΩG(cG1PGi,t+cG0)ΔtCDG=∑t∈ΩT∑i∈ΩWcWcur(PW,Prei,t−PWi,t)Δt+∑t∈ΩT∑i∈ΩPVcPVcur(PPV,Prei,t−PPVi,t)Δt

(2)

式中：CM为配电网向主网购电产生的费用；CLoss为配电网线路潮流产生的网损费用；CBat为储能装置在调度期间充放电所需要的运行费用；CG为微型燃气轮机的发电成本；CDG为风力发电（windturbinegeneration，WTG）和光伏发电（photovoltaicgeneration，PVG）弃风弃光的惩罚费用；CM+CLoss为配电网侧产生成本，CBat+CG+CDG为可控资源侧成本；ΩT为所有调度时段组成的集合；ΩM为配电网中与主网交换功率的耦合节点集合；Ωline为配电网拓扑中所有线路组成的集合；ΩC为配电网拓扑中接入储能装置的节点集合；ΩG为配电网拓扑中接入传统分布式发电机的节点集合；ΩW为配电网拓扑中接入WTG的节点集合；ΩPV为配电网拓扑中接入PVG的节点集合；t为调度时段；Δt为调度步长；cMt表示分时电价；cLoss为单位线损产生的费用；cdis为储能装置释放单位电量需要的费用；cch为储能装置充入单位电量需要的费用；cG1和cG0为燃气轮机发电成本函数的一次项系数和常数项；cWcur和cPVcur分别表示WTG和PVG削减单位发电量的惩罚系数；PMi,t表示节点i从主网购入的功率；Iij,t表示t时段下支路ij的电流；rij表示支路ij的电阻；Pchi,t和Pdisi,t分别表示节点i所接入的储能装置在t时段下的充放电功率；PGi,t表示所接入的燃气轮机在t时段下的发电功率；PW,Prei,t和PPV,Prei,t分别为WTG和PVG在t时段的发电功率；PWi,t和PPVi,t分别为WTG和PVG在t时段的实际流入节点i的功率。

1.2约束条件

1.2.1潮流平衡约束

传统潮流计算方法包括寻求局部最优解、近似线性化以及对潮流约束进行凸松弛等［17-18］等。其中，凸松弛技术可对潮流方程中的二次项和负二次项等非凸项元素进行处理，能够保证算法的高效性和寻优结果的最优性［18］。基于配网调度优化对潮流计算的准确性要求，本文采用二阶锥松弛（second-orderconeprogramming，SOCP）建立潮流平衡约束，并且本文将在仿真分析中对二阶锥进行误差研究，验证二阶锥松弛的有效性。对辐射型配电网，可列写其DistFlow形式的潮流方程，具体表示如下［17］：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪∑i∈Ωu(Pij,t−I2ij,trij)+∑Pj,t=∑z∈ΩvPjz,t∀j∈ΩBUS,∀t∈ΩT∑i∈Ωu(Qij,t−I2ij,txij)+∑Qj,t=∑z∈ΩvQjz,t∀j∈ΩBUS,∀t∈ΩT

(3)

U2j,t=U2i,t−2(rijPij,t+xijQij,t)+(r2ij+x2ij)I2ij,t,∀ij∈Ωjline,∀t∈ΩT

(4)

I2ij,t=P2ij,t+Q2ij,tU2i,t,∀ij∈Ωline,∀t∈ΩT

(5)

∑Pj,t=−PLj,t+PMj,t+PWj,t+PPVj,t+PGj,t−Pchj,t+Pdisj,t,∀j∈ΩBUS,∀t∈ΩT

(6)

∑Qj,t=−QLj,t+QMj,t+QWj,t+QPVj,t+QSVGj,t,∀j∈ΩBUS,∀t∈ΩT

(7)

式中：ΩuΩu为配电系统中以jj为末端节点的支路所对应的首端节点集合；ΩvΩv为配电系统中以jj为末端节点的支路所对应的末端节点集合；ΩBUSΩBUS为配电网络拓扑中所有节点组成的集合；ΩjlineΩlinej表示与节点jj衔接的支路集合；Pij,tPij,t和Pjz,tPjz,t分别为与节点jj相连的支路ijij和jzjz在tt时段的有功潮流；Qij,tQij,t和Qjz,tQjz,t分别为与节点jj相连的支路ijij和jzjz在tt时段的无功潮流；Uj,tUj,t和Ui,tUi,t为节点j和ii在tt时段的电压幅值；xijxij为支路ijij的线路电抗；Iij,tIij,t为支路ijij在tt时段的电流幅值；∑Pj,t∑Pj,t表示节点jj在tt时段的有功功率净注入量；∑Qj,t∑Qj,t表示节点jj在tt时段的无功功率净注入量；PLj,tPj,tL、PMj,tPj,tM、PWj,tPj,tW、PPVj,tPj,tPV、PGj,tPj,tG、Pchj,tPj,tch和Pdisj,tPj,tdis分别表示节点jj在tt时段的负荷有功功率需求、主网提供的有功功率、WTG提供的有功功率、PVG提供的有功功率、微型燃气轮机提供的有功功率、储能装置充电功率和储能装置放电功率；QLj,tQj,tL、QMj,tQj,tM、QWj,tQj,tW、QPVj,tQj,tPV和QSVGj,tQj,tSVG分别表示节点jj在tt时段的负荷无功功率、主网提供的无功功率、WTG提供的无功功率、PVG提供的无功功率和静止无功发生器（staticvargenerator，SVG）提供的无功功率。

式（3）~（5）中存在Uj,tUj,t、Ui,tUi,t、Iij,tIij,t的二次项和Ui,tUi,t的负二次项，若代入模型计算会出现非凸问题。可通过二阶锥处理将模型转化为凸优化问题，具体计算如下：

⎧⎩⎨⎪⎪Ut2,i=U2i,t,∀i∈ΩBUS,∀t∈ΩTIt2,ij=I2ij,t=P2ij+Q2ijUt2,i,∀ij∈Ωline,∀t∈ΩT

(8)

∥∥[2Pij,t2Qij,tIt2,ij−Ut2,i]T∥∥2≤It2,ij+Ut2,i,∀ij∈Ωline,∀t∈ΩT

(9)

式中：Ut2,iU2,it和It2,ijI2,ijt分别替代Ui,tUi,t和Iij,tIij,t的二次项。将式（8）和式（9）代入式（3）~（5）可得到下式：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪∑i∈Ωu(Pij,t−It2,ijrij)+∑Pj,t=∑z∈ΩvPjz,t∀j∈ΩBUS,∀t∈ΩT∑i∈Ωu(Qij,t−It2,ijxij)+∑Qj,t=∑z∈ΩvQjz,t∀j∈ΩBUS,∀t∈ΩTUt2,j=Ut2,i−2(rijPij,t+xijQij,t)+(r2ij+x2ij)It2,ij,∀ij∈Ωline,∀t∈ΩT∥∥[2Pij,t2Qij,tIt2,ij−Ut2,i]T∥∥2≤It2,ij+Ut2,i∀ij∈Ωline,∀t∈ΩT

(10)

1.2.2配电网安全性约束

在配电网线路及元件无故障的情况下，配电网安全性约束主要包括线路静态电流限制和静态电压限制［19］，具体约束见下式：

0≤It2,ij≤(Imaxij)2,∀ij∈Ωline,∀t∈ΩT

(11)

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪(Umin)2≤Ut2,i≤(Umax)2∀i∈ΩBUS,∀t∈ΩTUt2,i=(UBASE)2∀i∈ΩM,∀t∈ΩT

(12)

式中：ImaxijIijmax表示线路电流幅值最大允许值；UmaxUmax和UminUmin分别为节点电压上、下限；UBASEUBASE为参考节点电压值，当节点ii为与主网相连的节点时，其电压为参考电压。

1.2.3主网交换功率约束

本文的日前调度模型负荷规模相比分布式电源规模大，配电网电无富余的电能出售给主网，因此本文不考虑配网售电给主网获取收益的形式。配电网和主网之间的交互功率应满足以下约束：

0≤PMj,t≤PMmax,∀j∈ΩM,∀t∈ΩT

(13)

式中：PMmaxPmaxM为配电网与主网交换功率上限值。

1.2.4可控分布式电源出力约束

配电网中的可控分布式电源可在配电网电力供应不足时当作备用电源，在主网电价过高和负荷高峰期时能够减少配电网运行成本。可控分布式电源主要有燃料电池、微型燃气轮机等。本文考虑配电网中包含微型燃气轮机的情况，其运行约束条件应包括输出功率约束和爬坡约束［20］，具体如下所示：

0≤PGj,t≤PGmax,∀j∈ΩG,∀t∈ΩT

(14)

−PGDn≤PGj,t−PGj,t−1≤PGUp,∀j∈ΩG,∀t∈ΩT

(15)

式中：PGmaxPmaxG为微型燃气轮机的最大功率输出；PGDnPDnG为微型燃气轮机的滑坡功率变化限度；PGUpPUpG为微型燃气轮机的爬坡功率变化限度。

1.2.5储能装置运行约束

1）充放电功率约束

储能装置（energystoragesystem，ESS）在同一时段时只能处于充电状态或放电状态，且其充放电功率不能超出设定值，具体见如下约束：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪xCj,t∈{0,1}0≤Pdisj,t≤xCj,tPdismax0≤Pchj,t≤(1−xCj,t)Pchmax∀j∈ΩC,∀t∈ΩT

(16)

式中：xCj,txj,tC为ESS的状态变量，当xCj,txj,tC取值为1时，ESS处于放电状态，当xCj,txj,tC取值为0时，ESS处于充电状态，xCj,txj,tC变量保证ESS不能同时处于充放电的状态；PchmaxPmaxch和PdismaxPmaxdis分别为ESS交流侧的最大充、放电功率，其数值主要由储能并网逆变装置决定。

2）EES充放电容量约束

为有效延长ESS的使用寿命，调度决策方案中应禁止ESS过度充电或者过度放电。因此ESS在各个时段的电量应介于设备允许的上、下限之间［21］，具体见下式：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪ECj,t=ECj,0+ω∑t=1TPchj,tΔt−1ω∑t=1TPdisj,tΔtECmin≤ECj,t≤ECmax∀j∈ΩC,∀t∈ΩT

(17)

式中：ECj,tEj,tC为ESS在调度tt时段的剩余电量；ECj,0Ej,0C为调度阶段中ESS在初始时刻的电量；ECminEminC和ECmaxEmaxC分别为ESS在调度时允许的最小和最大容量；ωω为充放电效率系数；T为调度时段总数。

1.2.6DG相关运行约束

该部分约束主要包括WTG和PVG的弃风、弃光约束和无功功率控制约束［8］等，具体约束见下式：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪(1−εWmax)PW,Prej,t≤PWj,t≤PW,Prej,t(1−εPVmax)PPV,Prej,t≤PPVj,t≤PPV,Prej,t−PWj,tγtanθWLim≤QWj,t≤PWj,tγtanθWLim−PPVj,tγtanθPVLim≤QPVj,t≤PPVj,tγtanθPVLim∀j∈ΩW或ΩPV,∀t∈ΩT

(18)

式中：εWmaxεmaxW为风机有功出力的最大削减比例；εPVmaxεmaxPV为光伏有功出力的最大削减比例；PW,Prej,tPj,tW,Pre为风机在tt时段的有功出力；PPV,Prej,tPj,tPV,Pre为光伏在tt时段的有功出力；θWLimθLimW和θPVLimθLimPV分别为风机和光伏极限运行状态下的功率因数角，本文DG的额定功率因素均取0.85；为了避免DG在极限状态下运行，其无功出力需要设置阈值，γγ为阈值系数，本文取0.8。

1.2.7SVG运行约束

SVG运行约束具有快速、动态和连续调节无功的性能，能够及时提供配电网调度需要的无功功率，因此本文主要考虑SVG作为无功补偿装置。

SVG运行主要受其调控范围约束，具体见下式：

−QSVGmax≤QSVGj,t≤QSVGmax,∀j∈ΩSVG,∀t∈ΩT

(19)

式中：QSVGmaxQmaxSVG为SVG允许的功率输出最大绝对值。

1.2.8电价相关约束

在自由化的电力市场中电价由供需双方决定，并且由于电力商品具有不可贮存、实时平衡、输送约束等特性使得电价存在波动性［15］。电价向下波动时对优化结果的鲁棒性无影响［16］，故本文只考虑电价向上波动。

电价波动范围具体如下：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪cM,mint≤cMt≤cM,maxtcM,mint=cM,PretcM,maxt=(1+εM)cM,Pret∀t∈ΩT

(20)

式中：cM,PretctM,Pre为tt时段的预测电价；εMεM为电价预测误差因子；cM,mintctM,min为tt时段电价波动范围内的最小值；cM,mintctM,min为tt时段电价波动范围内的最大值。

式（2）中的CMCM包含双线性项cMtPMi,tctMPi,tM，使得优化模型无法直接求解。本文引入McCormick凸包络线法对其进行凸松弛［22］。需要说明的是，该模型的双线性变量只存在于目标函数，引入McCormick凸包络线法不会改变模型的可行域，从而保障优化模型不会因为松弛误差而求出违背运行约束的解集。此外McCormick凸包络松弛保证CMCM与cMtctM和PMi,tPi,tM严格遵守正相关关系，松弛后的误差也可通过仿真分析结果进行检验。松弛方法具体见下式：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪PM,maxi,t=σmax∑i∈ΩBUSPL,Prei,tPM,mini,t=σmin∑i∈ΩBUSPL,Prei,tPM,mini,t≤PMi,t≤PM,maxi,ti∈ΩM,∀t∈ΩT

(21)

ςi,t=cMtPMi,t

(22)

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪ςi,t≥cM,mintPMi,t+cMtPM,mini,t−cM,mintPM,mini,tςi,t≥cM,maxtPMi,t+cMtPM,maxi,t−cM,maxtPM,maxi,tςi,t≤cM,maxtPMi,t+cMtPM,mini,t−cM,maxtPM,mini,tςi,t≤cM,mintPMi,t+cMtPM,maxi,t−cM,mintPM,maxi,ti∈ΩM,∀t∈ΩT

(23)

式（21）~（23）中：PM,maxi,tPi,tM,max和PM,mini,tPi,tM,min分别表示tt时段节点ii与主网可能出现的最大交换功率和最小交换功率，两者相比tt时段的配电网有功负荷需求总量留有一定的裕度；σmaxσmax和σminσmin为负荷裕度系数，综合考虑负荷、分布式电源和储能装置的规模，本文σmaxσmax取3，σminσmin取0；ςi,tςi,t为McCormick凸包络方法引入的凸松弛变量，以此消除目标函数的双线性项。

2二阶段鲁棒日前调度模型及其求解方法

如果不确定性因素不存在波动性，配电网日前调度便是确定性规划问题，普通数学优化软件便可直接求解。

然而实际的配电网面临诸多随机因素的影响，预测精度难以保障，确定性调度模型将难以保证模型的解在不确定因素为“最恶劣”时仍为可行解，调度方案的准确性则取决于不确定因素的预测精度，使得调度方案具备风险性。

因此本文针对负荷需求、光伏出力、风电出力和实时电价的不确定性，提出了配电网的二阶段鲁棒日前优化调度模型，采用Benders分解法进行求解。

2.1二阶段鲁棒日前调度模型

2.1.1不确定集合

负荷需求、光伏出力、风电出力和实时电价的随机性需要考虑不确定参数的概率分布，但是其分布函数在现实中难以准确描述。盒式集合能够对不确定集合进行简化，并且保证模型的凸优化性质不会改变［23］。传统鲁棒方法的不确定变量仅取不确定区间的边界值，而盒式集合能够扩大可行域，具体表示为

ΩD=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪PˆLj,t,PˆLj,t∈[(1−τ)PL,Prej,t,(1+τ)PL,Prej,t],∀j∈ΩBUS,∀t∈ΩTQˆLj,t,QˆLj,t∈[(1−τ)QL,Prej,t,(1+τ)QL,Prej,t],∀j∈ΩBUS,∀t∈ΩTPˆW,Prej,t,PˆW,Prej,t∈[(1−τ)PW,Pre0j,t,(1+τ)PW,Pre0j,t],∀j∈ΩW,∀t∈ΩTPˆPV,Prej,t,PˆPV,Prej,t∈[(1−τ)PPV,Pre0j,t,(1+τ)PPV,Pre0j,t],∀j∈ΩPV,∀t∈ΩTcˆMt,cˆMt∈[cM,Pret,(1+εM)cM,Pret],∀j∈ΩM,∀t∈ΩT

(24)

式中：ΩDΩD为不确定因素的盒式集合；PL,Prej,tPj,tL,Pre和QL,Prej,tQj,tL,Pre分别表示节点jj在tt时段的有功和无功负荷的预测量；PW,Pre0j,tPj,tW,Pre0为节点jj在tt时段的WTG预期发电量；PPV,Pre0j,tPj,tPV,Pre0为节点jj在tt时段的PVG预期发电量；PW,Prej,tPj,tW,Pre和PPV,Prej,tPj,tPV,Pre分别为节点jj在tt时段WTG和PVG的实际发电量；cM,PretctM,Pre为主网在tt时段的预期售电价格，根据式（20）可知，电价波动只考虑向上波动；ττ和εMεM为鲁棒模型的不确定度。

不确定度可以表征鲁棒模型的随机性，其值越大，盒式集合出现的场景集合更“恶劣”，优化模型的随机性也更强。综上分析可知，盒式集合可由随机变量的基准值和不确定度确定。为简化分析，本文假设负荷需求、光伏出力、风电出力的不确定度都取ττ，电价的不确定度取εMεM。

2.1.2鲁棒优化数学模型

考虑鲁棒变量的不确定性，式（1）、（2）和（6）~（23）可以整理为以下形式：

F=minx∈{0,1},y(bTy)

(25)

Hy(d)≤h

(26)

Ax+By(d)≤g

(27)

Iyy(d)+Idd≤p

(28)

∥My(d)∥2≤Ey(d)

(29)

式中：xx为一阶段决策变量的列向量，其向量元素均为二进制，本文中其由储能装置的状态决策变量组成；dd为二阶段鲁棒变量组成的列向量，由式（24）内的鲁棒变量组成；除去以上变量，原模型剩余的变量组成向量yy；bTbT为目标函数对应的系数矩阵；HH和hh分别是整理式（11）~（15）、式（17）和（19）得到的系数矩阵和常数项向量，分别对应配电网安全性约束、主网交换功率约束、可控分布式电源出力约束、EES充放电容量约束、SVG运行约束；AA、BB和gg为整理式（16）得到的系数矩阵和常数项向量，对应的是充放电功率约束；IyIy、IdId和pp为整理式（10）、（18）和（23）得到的系数矩阵和常数项向量，对应的是潮流平衡约束、DG相关运行约束、电价相关约束；MM和EE由整理式（9）得到的二阶锥系数矩阵，对应的是二阶锥约束。二阶段向量可调整自身数据适应二阶段鲁棒变量dd的变化，参考文献［12］的方法，式（25）~（29）可进一步整理为以下二阶段形式：

minx{maxd∈ΩD[miny∈Ω(x,d)(bTy)]},s.t.x∈{0,1}

(30)

式中：Ω(x,d)Ωx,d为yy的可行域，即Ω(x,d)={Hy≤h,Ax+By≤g,Iyy+Idd≤p,∥My∥2≤Ey}Ωx,d=Hy≤h,Ax+By≤g,Iyy+Idd≤p,My2≤Ey。

第2阶段中存在max-min形式的优化问题，用目前优化领域的算法和优化求解器均无法直接进行求解，本文利用对偶理论和二阶锥对偶理论［13］对式（30）进行处理，具体形式如下：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪R(x)=maxφ,λ,η,u,v,d[λT(Ax−g)−φTh+ηTIdd−ηTp]φTH−λTB−ηTId+uTM+vTE=bTφ≥0,λ≥0,η≥0,v≥0∥u∥2≤vd∈ΩD

(31)

式中：R(x)Rx为第2阶段优化模型将min-max目标函数对偶后得到的新形式，其中xx在第2阶段的优化过程中视为已知变量；φφ、λλ和ηη为式（26）~（28）对应的对偶变量，uu和vv为式（29）对应的二阶锥对偶变量。

2.2鲁棒优化数学模型的求解方法

根据2.1.2节可知，二阶段鲁棒优化模型的第1阶段的目标是寻找最优的储能装置运行状态集，第2阶段的目标是在储能装置运行状态已知和不确定变量为“最恶劣”的情况下寻找最优的调度方案。

二阶段鲁棒模型可采用Benders分解法，将原问题分为主问题和子问题再进行迭代求解。主问题结合历代子问题获取的最优割平面进行求解，从而获取新的决策变量xx，接着将新的决策变量xx代入子问题进行求解，依次反复进行迭代，直到主‐子问题的目标函数达到收敛的条件，最终可得到原问题的最优解［12］。

子问题的目标函数中存在双线性项，本文使用OA法［24］进行求解。

2.2.1Benders分解法的求解过程

1）参数初始化

随机取满足第1阶段的约束条件初始解x0x0，将x0x0代入式（31）定义的第2阶段优化模型R(x)Rx，即子问题。求解R(x0)Rx0，并获取子问题的初始解(φ1,λ1,η1,u1,v1,d1)φ1,λ1,η1,u1,v1,d1。设置Benders主问题的上界值UBD=+∞UBD=+∞和下界值LBD=−∞LBD=-∞，设置当前迭代次数k=1k=1，设置Benders迭代过程中的允许误差参数ε(ε0)εε0。

2）求解Benders分解法的第kk代主问题

利用混合整数线性规划（mixedintegerlinearprogramming，MILP）对主问题进行求解，主问题形式如下：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪F=minx(α)α≥λlT(Ax−g)−φlTh+ηlTIddl−ηlTp,∀l≤kx∈{0,1}

(32)

式中：αα为Benders分解法的辅助变量；φlφl、λlλl、ηlηl和dldl为迭代次数为ll时子问题R(x)Rx对应的最优解。

根据式（32），求解主问题可得最优解(xk,αk)xk,αk，将主问题求得的目标函数值赋给Benders分解法的下界值，即LBD=αkLBD=αk。

3）求解Benders分解法的第kk代子问题

求解子问题R(xk)Rxk得到最优解集(φk+1,λk+1,ηk+1,uk+1,vk+1,dk+1)φk+1,λk+1,ηk+1,uk+1,vk+1,dk+1，将子问题求得的目标函数值赋给上界值，即UBD=R(xk)UBD=Rxk。

4）检验迭代终止的判断条件(UBD−LBD)≤εUBD-LBD≤ε是否成立

若判断条件成立，则迭代达到终止条件，此时xkxk鲁棒优化模型第1阶段的最优解，(φk+1,λk+1,ηk+1,uk+1,vk+1,dk+1)φk+1,λk+1,ηk+1,uk+1,vk+1,dk+1为第2阶段的最优解。若不成立，令k=k+1k=k+1，返回步骤（2），继续进行迭代。

2.2.2Benders分解法子问题的求解过程

式（31）中存在双线性项ηTIddηTIdd，使得优化模型无法直接求解。由于无法提前获知对偶变量ηη的取值范围，因此不能采用McCormick凸包络松弛。本文将利用OA法将问题分为OA主问题和OA子问题。OA法的求解步骤如下所述。

1）参数初始化

获取主问题的解集xkxk，随机选取某一场景的不确定变量d1∈ΩDd1∈ΩD，设置OA迭代过程的上界值UOA=+∞UOA=+∞和下界值LOA=−∞LOA=-∞，当前迭代次数n=1n=1，设置OA法迭代过程中的允许误差参数δ(δ0)δδ0。

2）求解OA子问题

第nn代OA子问题的具体形式如下：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪S(xk,dn)=maxφ,λ,η,u,v[λT(Axk−g)−φTh+ηTIddn−ηTp]φTH−λTB−ηTId+uTM+vTE=bTφ≥0,λ≥0,η≥0,v≥0∥u∥2≤v

(33)

在已知(xk,dn)xk,dn的情况下，对S(xk,dn)Sxk,dn进行求解，进而获取原调度优化对偶模型的最优解(φn,λn,ηn,un,vn)φn,λn,ηn,un,vn。将OA子问题求得的目标函数值赋给OA法的下界值，即LOA=S(xk,dn)LOA=Sxk,dn。定义OA法在第nn次迭代产生的最优割平面Ln(d,η)Lnd,η，其具体形式如下所示：

Ln(d,η)=ηnTIddn+(η−ηn)TIddn+(Idd−Iddn)Tηn

(34)

3）检验迭代终止的条件

判断(UOA−LOA)≤δUOA-LOA≤δ是否成立，若判断条件成立，则迭代达到终止条件，当前迭代的最优解(φn,λn,ηn,un,vn,dn)φn,λn,ηn,un,vn,dn是OA法求解原调度优化对偶模型获取的最终解。若不成立，令n=n+1n=n+1，进入OA法的步骤（4）。

4）求解OA主问题

第nn代OA主问题的具体形式如下：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪U(dn−1,ηn−1)=maxφ,λ,η,u,v,d,β[λT(Ax−g)−φTh+β−ηTp]β≤Li(d,η),∀i≤n−1φTH−λTB−ηTId+uTM+vTE=bTφ≥0,λ≥0,η≥0,v≥0d∈ΩD,∥u∥2≤v

(35)

将历代获取的最优割平面Li(d,η)Lid,η代入OA主问题，此时U(dn−1,ηn−1)Udn-1,ηn-1是线性优化问题，利用Cplex数学优化软件便可进行求解。然后将OA主问题求得的目标函数值赋给上界值，即UOA=U(dn−1,ηn−1)UOA=Udn-1,ηn-1。将求得的dndn代入步骤（2），继续进行迭代。

具体的鲁棒优化数学模型的求解逻辑详见图2。

图2鲁棒优化数学模型的求解逻辑图

Fig.2Logicdiagramforrobustoptimizationmathematicalmodelsolution

3算例仿真分析

3.1算例说明

本文以改进的电气和电子工程师协会（InstituteofElectricalandElectronicsEngineers，IEEE）33节点配电系统［25］为例进行仿真分析，设调度周期为24h，单位调度时间间隔为1h；节点额定电压取12.66kV，各个节点的电压波动范围为额定电压的0.95~1.05倍；线路额定电流取0.5kA；在节点1处配电网与主网耦合，在节点3、6和11处接入微型燃气轮机，在节点18和22处接入WTG，在节点25和33处接入PVG，在节点2、17、21、24和32处接入储能装置；配电网的拓扑图详见图3。

图3算例拓扑结构

Fig.3Thetopologyoftheexample

SVG的可调范围为-500~500kVar；配电网与主网交互功率限制为3000kW；算例的网损成本取0.68元/kWh；主网分时电价参考文献［11］，具体见图4；风光出力和配电网日负荷见参考文献［8］，为便于数据分析，本文假设每个节点的负荷时序一致，曲线详见图5，图中每台WTG有功出力的基准值为100kW，每台PVG有功出力的基准值为100kW，配网总负荷有功功率基准值为1980kW，总负荷无功功率基准值为990kVar；求解算法中，Benders分解法迭代最大相对误差设为0.5%；OA法迭代最大相对误差设为0.2%；盒式集合中的不确定度ττ取20%，εMεM取20%。WTG［7］的成本系数cWcurccurW和PVG［7］的成本系数cPVcurccurPV均为0.544元/（kW⋅h）；储能装置的运行参数［11］PdismaxPmaxdis和PchmaxPmaxch均为500kW，EC0E0C、ECminEminC、ECmaxEmaxC分别为1000、400、1800kW⋅h，充放电成本系数cdiscdis和cchcch均为2.47元/（kW⋅h），充放电效率系数ωω为0.95；微型燃气轮机［11］的最大输出功率PGmaxPmaxG为300kW，成本系数cG1c1G和cG0c0G分别为4.36元/（kW⋅h）和0元/（kW⋅h），爬坡约束功率PGUpPUpG和PGDnPDnG均为200kW。

图4分时电价

Fig.4Time-sharingprice

图5负荷需求、PVG出力和WTG出力的预测曲线

Fig.5Predictioncurvesofloaddemand,PVGoutputandWTGoutput

3.2仿真结果及其分析

3.2.1调度优化方案的合理性分析

调度优化结果如图6~8所示。图6为各台微型燃气轮机发电功率和配电网向主网购入的有功功率分布情况。图7为各台储能装置在各个时段的充放电功率分布情况，充电时对应柱状图颜色为绿色，放电时对应柱状图为蓝色。图8为SVG各个时段提供的无功功率分布情况，左侧坐标轴对应的是折线图，右侧纵坐标对应的是柱状图。

图6燃气轮机输出功率及主网购电功率

Fig.6Outputpowerofgasturbineandpowerpurchasedfrommainnetwork

图7储能装置运行状态

Fig.7Operationstatusofenergystoragedevice

图8SVG的无功功率输出分布情况

Fig.8ReactivepoweroutputdistributionofSVG

由图4~6可知，微型燃气轮机只在主网电价高峰期（11~15h和19~21h）和平缓期（8~10h、16~18h和22~23h）时进行发电。由此可见，本模型能有效安排微型燃气轮机的功率分布，减少配电网的运行成本。由图7可知，所有储能装置能够在主网电价低峰期时及时充电，而后在电价高峰期提供电力，有效协助配电网“削峰填谷”。为便于数据分析，本文中的仿真模型已经假设各个节点的负荷时序特性一致。由图8可知，各个SVG无功出力的时序特性与负荷的无功时序特性基本一致，且配电网中大部分无功功率由SVG提供。这说明SVG能够就近为负荷提供无功功率，减少无功功率输送的距离，因此SVG能够达到支撑电压降低和减少线路损耗的作用。

3.2.2鲁棒优化模型的有效性分析

为说明鲁棒优化模型的有效性，本文同时对不考虑鲁棒性的确定性模型和不同不确定度的模型进行求解，其日前调度成本详见表1。

表1配电网日前调度成本

Table1Costofday-aheaddispatchingofdistributionnetwork

优化模型ττεMεM配电网侧费用/元可控资源侧费用/元总费用/元主网购电

费用

网络损耗

费用

储能运行

费用

燃气轮机

发电费用

不考虑鲁棒性的确定性模型--37401.22145.5013538.5759911.67110996.97不同不确定度的模型000.050.0538497.24145.4313538.7859911.68112093.15000.20.239680.60140.3113539.0763813.89117173.890.050.050041967.22170.0416217.3560071.35118425.970.20.20069385.70267.1218509.1960168.00148330.030.20.20.20.276002.35258.8618509.2065400.00160170.42

根据表1可知，随着预测不确定度的增加，配电网的运行成本也不断提高。这说明鲁棒经济调度模型会根据调度者的冒险意愿，牺牲一定的经济性来保证配电网对更恶劣场景的适应性。为进一步说明本模型鲁棒调控的有效性，本文进行以下分析。

1）调度成本与电价预测的不确定度εMεM之间的关系

随着电价预测的不确定度εMεM的增加，储能运行费用一直保持不变，这是因为在负荷需求不变的前提下，储能只须在电价低谷时存储恒定的电量，而后将这些电量在电价高峰时完全释放。主网购电费用因为电价整体抬高而保持上升状态。

燃气轮机发电费用先维持不变，再呈现上升状态。根据表1和图4可知，配电网在电价低谷期时从主网购电的成本最低，而在电价平缓期和高峰期时燃气轮机发电成本比主网电价更低，此时应保持高功率发电状态。然而燃气轮机的爬坡约束导致其无法迅速改变发电功率，这使得燃气轮机是否在电价低谷时应提前进入发电状态以适应非电价低谷期时的高功率发电状态存在争议。燃气轮机在电价低谷时为下一调度时段的高功率发电状态所额外付出的最小费用称为“缓冲成本”。当燃气轮机的发电成本比主网购电成本低时，因爬坡约束无法供应足够发电量而被迫选择从主网购电进行替代而产生的费用称为“不可控效益”。随着电价整体不断抬高，“不可控效益”随之上涨，“缓冲成本”则保持不变。只有当“不可控效益”比“缓冲成本”高时，燃气轮机才会在电价低谷期提前进入发电状态，使得燃气轮机发电费用额外包含“缓冲成本”。

2）调度成本与不确定度ττ之间的关系

随着不确定度ττ的增加，运行成本中的各项费用都随之上涨。这是因为配电网为了满足更多的负荷需求，各时段的购电量和微型燃气轮机的发电量均呈上升趋势。因为容量有限，储能装置在电价低谷期电量贮满时便无法继续充电，因此其运行费用在不确定度ττ增加到一定数值时便不再随之升高。

3.2.3配电网不同配置对鲁棒优化结果的影响

由上文分析可知，合理安排微型燃气轮机和储能装置的功率分布能有效减少配网运行成本。为了给配电网规划者提供指导性意见，本文将分析在不同不确定度下微型燃气轮机的最大输出功率和储能装置容量与调度成本的灵敏度关系，具体如图9和图10所示。图9（a）和图10（a）的电价不确定度εMεM为0.2，图9（b）和图10（b）的不确定度ττ为0.2。

图9燃气轮机最大输出功率与调度成本的灵敏度关系

Fig.9Sensitivityrelationshipbetweenmaximumoutputpowerofgasturbineandschedulingcost

图10储能装置容量与调度成本的灵敏度关系

Fig.10Sensitivityrelationshipbetweencapacityofenergystoragedeviceandschedulingcost

分析图9~10的曲线变化趋势可知，随着微型燃气轮机的最大功率和储能装置容量的增加，配电网的调度成本先逐渐减少，后趋于平缓发展，并且不确定度越高，配电网的调度成本下降越快。这说明选择性能越好的配电网配置能更好地调节模型鲁棒性，且不确