考虑风电相关性的电力系统随机无功备用优化

方斯顿，程浩忠，宋 越，曾平良，姚良忠

（1.上海交通大学 电力传输与功率变换控制教育部重点实验室，上海 200240；2.中国电力科学研究院，北京 100192）

0 引言

实践表明，无功备用对维护电力系统稳定运行有着重要价值[1-5]。对系统的无功备用水平进行分析进而提出优化方法是改善运行安全水平的有效手段。目前，稳态情况下的无功备用优化评估方法已有不少成果，主要可分为物理备用[6-7]和有效备用[8]，两者的区别在于前者将发电机的物理无功上限与基态无功出力之差作为备用值，而后者则将无功备用定义为系统临界点处的无功出力与基态出力之差。研究表明，在系统临界点处大多数发电机的无功出力均没有越限[9]，采用物理备用会乐观地估计系统的无功水平，因此采用有效备用来衡量发电机的无功支撑能力更加准确。而由于发电机的安装位置不同，不同发电机的无功备用价值不尽相同，将每台发电机的无功备用简单相加并不能准确衡量系统的总无功备用。因此有必要将系统的总无功备用定义为各台发电机无功备用的加权和。文献[10]以发电机无功出力对潮流雅可比矩阵最小奇异值的灵敏度定义权系数，但当系统规模扩大时，发电机无功出力对最小灵敏度的影响很小，因此该方法不适合大系统分析。文献[7-8]采用节点无功参与因子定义无功备用权系数，但仅适用于鞍结分岔SNB（Saddle Node Bifurcation）类型而不适合于极限诱导分岔LIB（Limit Induced Bifurcation）类型。文献[11]指出极限曲面法向量能够指示LIB时节点电压的薄弱性，但由于极限曲面的复杂性，法向量指标在某些区域可能变化较剧烈，不适合单独使用。

随着以风能为代表的新能源持续接入，节点注入功率的不确定性日益明显。目前在无功备用优化领域，考虑随机性的模型还未见报道。但众多研究表明，在含不确定性的系统中直接应用确定型模型的结果会增加系统的运行风险[12]。此外，文献[13]指出考虑风电相关性会恶化系统运行环境。因此将风电的随机性和相关性纳入无功备用优化模型是值得尝试的研究思路。

针对以上问题，本文定义了一种基于节点参与因子和极限曲面法向量乘积的无功备用权系数，并在模型中考虑风电随机性与相关性，提出电力系统随机无功备用优化模型，采用混合蛙跳算法SFLA（Shuffled Frog Leaping Algorithm）进行求解。对9节点系统和IEEE 118节点系统的仿真验证了本文方法的有效性。

1 模型的建立

1.1 发电机无功备用权系数

发电机无功备用权系数表征在系统不同位置的发电机对系统无功贡献的大小，贡献大的权系数也大，即在优化过程中价值较高，需要优先提高。针对目前权系数定义的缺陷，本文提出采用基于节点参与因子和极限曲面法向量指标乘积的权系数定义方式。

1.1.1 基于极限曲面法向量的节点电压薄弱性指标在节点注入功率空间中，极限运行曲面∑定义为：

其中，P、Q 分别为有功、无功节点注入功率；f（·）=0为临界点处的潮流方程；ΩSNB、ΩLIB分别为SNB和LIB电压稳定临界点集。根据文献[11]，∑的法向量在张成空间span（Q）上的投影nQ可以表征节点电压的薄弱性，nQ越大，代表该节点对无功的扰动越敏感，越需要增加无功支撑。

1.1.2 基于节点参与因子的节点电压薄弱性指标

参与因子法是一种特征值分析法[14]，设系统共有nb个节点，首先对潮流雅可比矩阵进行特征值分解，得到：

其中，ΔP、ΔQ、Δθ 和 ΔU 分别为节点注入有功、无功、节点电压相角和幅值的微分，均为nb-1维列向量；矩阵A和BT分别为左、右特征向量对应的矩阵。

由式（2）易得：

其中，Ai和Bi分别为A和B的第i列。

令 ΔP=0、ΔQ=ξj（即第 j个量为 1、其余为 0 的单位向量），可得第j个节点电压灵敏度为：

其中，Anb+j-1，i、Bnb+j-1，i分别为 A 和 B 中第 nb+j-1 行、第i列的元素；pj，i为状态j对模态i的参与因子。由式（4），系统中对应特征值最小的模态对系统的行为影响最大，而参与因子较大的节点对的贡献也大，也即该节点越需要无功支撑。不失一般性，下文的参与因子均针对主模态（即最小特征值对应的模态），将此模态下第j个节点的参与因子表示为pj。

1.1.3 本文定义的无功备用权系数

在nQ和pj，i的基础上，本文定义的发电机无功备用权系数，其元素为：

其中，wGi为发电机i的无功备用权系数；Ngen为系统的发电机集；NPQ为系统的PQ节点集；为节点j电压对发电机节点i端电压灵敏度，其值越大代表发电机主动控制节点电压的能力越强，为发电机节点i的无功出力对于节点j无功扰动的灵敏度，其值越大代表发电机响应节点无功扰动的能力越强，这2项衡量发电机对节点电压的控制能力和无功支撑能力；pj为节点j对主失稳模式的参与因子，nQj为nQ中节点j对应的元素，如前文所述，这2项衡量节点电压的薄弱性。

为避免数值差距影响指标准确性，本文采用将上述4项相乘的形式。这一指标综合表征了发电机对节点电压的控制能力及无功响应能力，其值越大代表其对薄弱节点的无功支撑能力越强，在无功备用优化中的价值越高。为便于分析比较，下文将不包含nQj项的权系数称为基于参与因子的权系数，而将不包含pj项的权系数称为基于法向量的权系数。

1.2 随机无功备用优化模型

a.目标函数如下：

其中，E（·）表示求期望值；QGR为本文定义的全系统总无功备用；wGT为1.1节中定义的无功备用权系数；QGc=[QGc1，…，QGci，…]T和 QG0=[QG01，…，QG0i，…]T分别为系统发电机在临界点和基态情况下的无功出力向量。

b.变量约束为：

其中，f为潮流方程；U和θ分别为节点电压幅值及其相角；Pl和Ql分别为节点的有功负荷以及无功负荷，风电机的有功、无功出力也作为负的负荷计入这一部分；QG0为发电机的基态无功出力；u为控制变量，包括发电机端电压、节点无功补偿量、可调变压器变比，为简化计算，本文控制变量均作连续化处理；下标max和min分别表示上限和下限；Pr（·）表示事件发生的概率；p为给定的越限概率。

2 模型的求解

前文所述的模型是带机会约束的非线性规划问题，本文采用随机潮流和混合蛙跳算法求解该问题。

2.1 基于均匀采样的随机潮流

均匀抽样法不同于普通的采样方法，其以确定性方式在采样空间中生成服从均匀分布的序列，可提高蒙特卡洛的精度[15]。过程简述如下，流程见图1。

b.获得m维独立满足多项式分布C（N，1/N）的向量 B1×m=[b1，b2，…，bm]。

c.利用简单随机采样获得在[-0.5，0.5]满足均匀分布的采样矩阵 Em×N=[e1，e2，…，em]T，其中 ek=[ek1，ek2，…，ekN]（k=1，2，…，m）。

d.采样矩阵 Ym×N=[y1，y2，…，ym]T可由式（11）生成，其中 yk=[yk1，yk2，…，ykN]（k=1，2，…，m）。

其中，i=1，2，…，m；j=1，2，…，N；Mod（·）表示取余运算。

图1 蛙跳算法求解无功备用优化流程Fig.1 Flowchart of reactive power reserve optimization by SFLA

获得Ym×N后，可根据文献[16]中的方法首先获得含相关性的正态分布样本，最后转化为所需的分布。最终采用蒙特卡洛获得输出变量的数字特征。

2.2 基于混合蛙跳的优化算法

混合蛙跳算法由Eusuff和Lansey在2003年提出[17]，在求解大规模优化问题时具有收敛速度快、全局寻优能力强的优点，在各学科中均得到应用[18-19]。现将其简要介绍如下。

以系统某一无功出力方案作为一只青蛙Xf=[xf1，xf2，…，xfnc]，其中 nc为系统中控制变量个数，控制变量包括静态电容器出力、发电机端电压、变压器分接头。

步骤1：设定初始参数。设定群体青蛙数量F、种群数mp、局部搜索迭代次数Nl、全局迭代次数Ng、青蛙所允许单次最大移动距离Smax。

步骤2：随机产生F只青蛙作为初始蛙群P=｛X1，X2，…，XF｝，第 f只青蛙代表问题的一个解 Xf=[xf1，xf2，…，xfnc]。

步骤3：对每只青蛙对应的方案进行连续潮流计算，获得QGc；按风速分布产生含相关性的风速序列，并对该方案进行随机潮流计算，获得该方案对应系统的无功备用期望值及节点电压、支路潮流和发电机无功出力越限的概率。适应度按式（12）计算。

其中，F0为给定正数以保证足够的选择压力；a、b、c均为大于1的常数，a=b＜c。

步骤4：若违反约束的青蛙数大于等于F/3，返回步骤2，重新生成青蛙替换不可行的青蛙，直到种群中2F/3以上的青蛙是可行的。

步骤5：划分青蛙种群。对所有青蛙按适应度降序排列，记适应度最好的青蛙为Pg，并将青蛙均分为mp个子蛙群，记每个蛙群数目为q。其中，第1只青蛙放入种群1，第mp只放入种群mp，第mp+1只放入种群1，依此类推。

步骤6：局部搜索。对于每个种群，记适应度最好的青蛙为 Pb，最差的为 Pw。 如式（13）、（14）对 Pw逐位进行更新得到P*w，若优于Pw，则将其替换；若未改进，使用Pg替换Pb，重复这一过程直到P*w优于Pw或达到局部搜索次数上限Nl；在更新过程中若有青蛙越限，则将其放置于边界上继续进行更新。

其中，pb、pw和 p*w分别为 Pb、Pw和 P*w的某一位。

步骤7：混合蛙群，重复步骤5、6，直到达到全局迭代次数Ng。

3 算例分析

3.1 9节点系统

为说明本文所提无功备用权系数的有效性，首先以9节点系统为例，详细数据参见Matpower 4.1[20]。系统接线图如图2所示，节点1为平衡节点，节点2、3为PV节点，其余为PQ节点。发电机基态情况与临界状态下系统潮流数据如表1所示，表中数据均为标幺值。临界状态下节点1发生无功越限，电压稳定裕度为1.3055，电压失稳类型为LIB。

图2 9节点系统单线图Fig.2 Single line diagram of 9-bus system

表1 基态和临界潮流解Table 1 Data for initial state and critical state

为比较本文定义的权系数、基于参与因子的权系数、基于法向量的权系数的效果，本文首先设计以下3种提高发电机无功备用的调度方案，各方案计算结果见表2，表中Δλ/λ表示电压稳定裕度增加百分比。

a.方案A：发电机1无功上限提高0.5 p.u.。

b.方案B：发电机2无功上限提高0.5 p.u.。

c.方案C：发电机3无功上限提高0.5 p.u.。

表2 3种指标与各方案下电压稳定裕度增量Table 2 Augment of voltage stability margin by 3 methods for 3 schemes

从表2中数据可知，本文定义的权系数能够更准确地区分不同发电机的无功备用价值。为进一步说明，在方案A和方案B同时实施的基础上本文再设计3种调度方案，计算结果见表3，表中数据均为标幺值。此时无功备用权系数分别为0.4798、0.1770、0.0754，调度前电压稳定裕度为1.3515，失稳类型为SNB。

a.方案D：节点1端电压调节至1.05 p.u.。

b.方案E：节点2端电压调节至1.05 p.u.。

c.方案F：节点3端电压调节至1.05 p.u.。

表3 各方案下发电机无功备用变化情况Table 3 Data of generator reactive power reserve for 3 schemes

3.2 IEEE 118节点系统

3.2.1 算例设置

为说明本文所提随机无功备用优化方法的有效性，以Matpower 4.1[20]中的IEEE 118节点系统为例对优化算法进行分析，该系统共有54个发电机节点、9台变压器和17个电容补偿节点，补偿容量上限均为100 Mvar。节点电压上限均设为1.06 p.u.，下限为0.96 p.u.，置信概率均设为 98%。在节点 50、51、57、58、93、94、95、96 这 8 个节点装设 200 MW 的风电场，其中 50、51、57、58 为同一区域，93、94、95、96 为另一区域，区域内风电相关性系数为0.8，区域间无相关性，风速均满足威布尔分布W（10.7，3.97）。系统的负荷分为3个区域：

a.区域 A：节点 24、33—73、116。

b.区域 B：节点 74—112、118。

c.区域 C：节点 1—23、25—32、113—115、117。

区域内的负荷相关系数均为0.5，不同区域内负荷相关系数为0。设负荷均满足正态分布，期望值等于负荷初始值，标准差为其期望值的5%。

算法参数设置为：青蛙数量F为200，局部搜索次数N1为10，全局搜索次数Ng为200，青蛙单次最大移动距离 Smax为 5，a=b=2，c=3，子蛙群数目 mp为10。均匀采样规模为500，待采样随机变量数为107（99个负荷+8个风电场）。

3.2.2 本文方法优化结果分析

优化前节点电压2%分位数和98%分位数如图3所示，图中节点电压为标幺值，后同。

图3 优化前节点电压分位数曲线Fig.3 Quantile curves of node voltage before optimization

从图3中可以看出，考虑风电随机性之后，优化前运行点违反机会约束14个，运行的风险大。优化后结果如图4所示，优化后节点电压的方差显著减小，机会约束均得到有效满足，系统的运行风险降低。为进一步说明本文方法相对确定型无功备用优化的优势，优化前后节点114电压的概率分布曲线如图5所示。

由图5可见，优化前节点114电压的期望值为0.9604 p.u.，满足节点电压约束下限，但考虑风电随机性后该节点电压越下限概率达到37%。因此忽略风电随机性会增加系统运行的风险。而采用本文方法使得节点114电压的方差减小，概率分布很好地满足了机会约束。

图4 优化后节点电压分位数曲线Fig.4 Quantile curves of node voltage after optimization

图5 优化前后节点114电压概率分布Fig.5 Voltage probability distribution of node 114 before and after optimization

在严重故障方面，经过N-1故障筛选，支路8-5开断是严重故障，电压稳定裕度仅0.67，远远低于基态情况。而本文方法得到方案在支路8-5开断后的负荷裕度由2821.79 MW上升到3651.86 MW。优化前、后系统节点电压2%分位数如图6所示。

图6 支路8-5开断时本文方案与基态电压2%分位数曲线Fig.6 2%quantile curve of node voltage for proposed and original methods when branch 8-5 is open

由图6可见，本文方法可以显著改善支路8-5切除后系统的电压水平，仅有个别节点电压违反机会约束。这说明本文方法不仅可以改善正常运行情况下的电压水平，降低系统运行风险，对故障情况下的电压水平和稳定裕度也有提高作用。

3.2.3 3种无功备用定义方式对比分析

为分析3种无功备用定义方式的有效性，3种无功备用定义方式下系统最薄弱节点95优化前、后的PV曲线如图7所示，图中λ为负荷增长因子。本文定义方式下电压稳定裕度由7989.37 MW上升到8657.78 MW，增加8.37%，远高于基于参与因子的1.06%和基于法向量的3.19%。这再次证明本文的无功备用定义方式可有效区分不同发电机的无功价值。

图7 优化前后节点95 PV曲线Fig.7 PV curve of node 95 before and after optimization

3.2.4 算法效率分析

为衡量算法寻优的效率及鲁棒性，在相同种群规模下，分别采用粒子群优化PSO（Particle Swarm Optimization）算法、遗传算法 GA（Genetic Algorithm）和混合蛙跳算法求解1.2节中模型100次，所得平均无功备用、最大无功备用QGRmax、最小无功备用QGRmin、达到QGRmax次数NQ、平均迭代次数N2、平均计算时间t列于表4中。

表4 各算法收敛特性Table 4 Convergence characteristics for different algorithms

为分析种群混合策略对混合蛙跳算法寻优的影响，表中SFLA1代表包含种群混合策略的混合蛙跳算法，而SFLA2不包括种群混合，种群数mp=1。遗传算法参数：迭代次数200，染色体数100，交叉率0.8，变异率0.1，采用轮盘赌选择机制。粒子群优化算法参数为：迭代次数200，粒子数100，惯性权重0.8，局部与全局学习因子2，粒子最大移动距离5。计算仿真平台基于 MATLAB2013a，Intel Core dual i7-3820 3.6 GHz，RAM 8 GB。

从表中可看出，4种算法所得结果相近，但SFLA1达到QGRmax次数NQ明显高于其他3种算法，平均迭代次数N2则低于其他算法，这证明本文采用的混合蛙跳算法具有较高的寻优效率。而SFLA2与粒子群优化算法所得结果相近，说明SFLA1能够高效寻优的原因在于其不仅采用类似粒子群的局部搜索策略，还通过混合不同种群个体，极大增加了全局寻优能力。此外，由于引入种群混合策略，SFLA1单次迭代的时间高于其他3种方法，但由于平均迭代次数N2小，最终计算时间仍然优于其他3种方法。

4 结论

本文提出了一种考虑风电相关性的电力系统随机无功备用优化方法，并定义了新的发电机无功备用权重因子。9节点算例的分析证明了本文无功备用定义的有效性，IEEE 118节点系统的仿真结果表明本文的优化算法有效，优化方案显著提高了系统的电压稳定裕度与电压水平，并降低了系统的运行风险。得到结论如下：

a.本文方法能够在SNB和LIB这2种情况下准确识别出系统中的关键发电机，识别效果优于基于参与因子的定义方式和基于法向量的定义方式，基于本文定义的无功备用优化模型能够有效提高系统的电压稳定裕度；

b.系统不同发电机的无功备用价值随安装位置不同变化很大，合理增大关键发电机的无功备用可以快速有效提高系统的电压稳定性裕度；

c.考虑风电相关性会使部分状态变量违反机会约束，增大系统的运行风险；

d.本文方法可以有效提高系统电压稳定裕度，改善系统的无功电压水平，减小由风电接入带来的电压波动，降低系统的运行风险，对改善严重故障下系统的电压稳定裕度有一定作用；

e.本文采用的混合蛙跳算法采用种群混合机制，寻优效率高，高于同样本规模的粒子群优化算法和遗传算法。

为简化计算，本文对电容补偿量、变压器分接头等离散变量均进行了连续化处理。此外在考虑负荷增长的场景时仅考虑等比例增长。以上问题尚有待进一步研究。