混合整数线性规划形式的抗差状态估计方法

陈艳波，马 进，陈 茜

（1.华北电力大学 新能源电力系统国家重点实验室，北京 102206；2.悉尼大学 电气与信息学院，澳大利亚 悉尼 NSW2006）

0 引言

电力系统状态估计SE（State Estimation）是能量管理系统EMS（Energy Management System）的基础和核心，是电网安全运行的重要保障[1-2]。SE算法最早由 Schweppe 和 Wildes在 1970 年提出[3]，至今已有多种SE方法。目前在国内外广泛应用的SE方法是加权最小二乘法 WLS（Weighted Least Squares）[3]及在此基础上提出的快速分解法FDSE（Fast Decoupled State Estimation）[4]。 WLS 和 FDSE 的优点是原理简单、计算效率高；缺点是抗粗差能力弱，易受不良数据的影响[5]。为此一般需要在WLS和FDSE之后加上一个基于残差的不良数据辨识环节，如最大标准化残差 LNR（Largest Normal Residual）方法[6-7]或估计辨识法[8-10]等。LNR方法和估计辨识法对单一不良数或弱相关的多不良数据具有较好的辨识能力，而对强相关的多不良数据（包括杠杆点不良数据）的辨识能力不强[6-7]。

除WLS和FDSE之外，人们还提出了抗差状态估计 RSE（Robust State Estimation），在估计过程中自动抑制不良数据。已提出的RSE方法主要包括加权最小绝对值估计 WLAV（Weighted Least Absolute Value）[11-12]、二次线性估计 QL（Quadratic-Linear）、二次常数估计 QC（Quadratic-Constant）[13-15]等。近年来提出的以合格率最大为目标的状态估计MNMR（Maximum Normal Measurement Rate）[16-18]、指数型目标函数状态估计MES（Maximum Exponential Square）[19]、最大指数绝对值状态估计 MEAV（Maximum Exponential Absolute Value）[20-21]、t型抗差状态估计[22]、电力系统线性状态追踪方法[23]、双线性抗差状态估计方法[24]以及量测噪声自适应抗差状态估计方法[25]均属于 RSE 的范畴。

最近国外学者提出了一种混合整数非线性规划MINP（Mixed Integer Nonlinear Programming）形式的RSE 方法[26]（以下简称为 MINP）。 MINP 的优点是对包括杠杆点不良数据在内的多不良数据具有很好的抑制能力。MINP的缺点是：在数学上，MINP需要求解非线性非凸优化问题，采用基于梯度的方法求解时，从理论上难以保证获得全局解；MINP的计算效率很低，难以满足大规模网络的实时计算要求。以上缺点可归因于SE量测方程的非线性，若能采用某种变换将SE的原始非线性量测方程精确变换为线性化的量测方程，则以上2个缺点有望得到解决。

由国外学者提出的双线性状态估计BSE（Bilinear State Estimation）方法[27-29]给出了 SE 量测方程的精确线性化方法。若将BSE中的精确线性化量测方程应用于MINP，则其变为混合整数线性规划 MILP（Mixed Integer Linear Programming）模型。MILP不仅具有良好的抗差性和收敛性，而且从数学上可保证获得全局最优解，并具有较高的计算效率。

1 MINP概述

一般地，SE的非线性量测方程可表示为：

其中，zєRm为量测向量，一般包括支路功率量测、节点注入功率量测及节点电压幅值量测等，m为量测量的总个数；xєRn为包括所有节点电压幅值和相角（参考节点相角除外）的状态变量向量，本文简称为状态向量，n=2N-1，N为网络中节点的总数目；为状态向量到量测向量的非线性映射，即量测表达式；εєRm为量测误差向量。

文献[26]基于量测不确定度理论，认为每个量测量均具有不确定性，任意量测量均可看作位于以量测值为中心的区间中，用公式可表示为：

其中，hi（x）为h（x）的第 i个分量；zi为 z的第 i个分量；ti-和ti+分别为用以衡量量测不确定度的下界和上界，它们可根据量测仪表的性能评定。

文献[26]认为，若所有的量测量均位于其容许的不确定区间内，则存在一个状态向量估计值，使得所有的量测量均满足式（2）；然而由于量测量中不可避免会含有不良数据，它们对应的量测估计值为状态向量估计值）不满足式（2），即此时不存在一个状态向量而使得所有的量测量均满足式（2）。为了将所有的量测量（包括正常量测和不良数据）表示为统一的不确定形式，文献[26]引入如下表达式：

其中，M为充分大的标量；bi为量测量zi对应的0／1二值变量，对于正常量测bi=0，对于不良数据bi=1。

基于式（3），文献[26]提出的 MINP模型如下：

MINP模型（4）的意义是：通过估计一个状态向量，使得尽可能多的量测量能够相互匹配。

2 MILP方法

2.1 MILP提出的动机和思路

MINP问题的根本原因在于量测方程的非线性（如式（1）所示），如果能采用某种变换将非线性量测方程（1）精确变化为线性量测方程，则MINP模型（4）将变为MILP问题，则全局寻优性和计算效率等问题均可得到解决。

文献[27-29]通过引进辅助状态向量和辅助量测向量，得到了如下形式的量测方程：

其中，y=｛Ui，Kij，Lij｝єRN+2b为辅助状态向量，Ui=Um2i为电压幅值的平方项，Umi为节点i的电压幅值，Kij=UmiUmjcosθij和 Lij=UmiUmjsinθij为支路 ij（i和 j为支路两端节点）对辅助状态向量的贡献，θij=θi-θj为支路 ij两端的相角差；为辅助量测向量，Pij和Qij分别为支路ij的有功和无功量测，Pi和Qi分别为节点i的注入有功和注入无功量测，Ii2j为支路ij的电流幅值量测的平方项，Iij为支路ij的电流幅值量测；BєRm×（N+2b）为常数雅可比矩阵，其非零元由支路（包括变压器）的导纳参数和变比决定；为噪声矢量；b为支路的总数目（并联支路等效为一条支路处理）。

BSE第1阶段基于式（5）构建了WLS模型，可得到辅助状态向量y的估计值，进一步利用其第2和第3阶段得到SE原始状态向量x的估计值。

对比式（1）和式（5）可发现，辅助状态向量和辅助量测向量的引入使得SE的原始非线性量测方程（1）变为精确线性化的量测方程（5），若将其应用于MINP，则可解决MINP存在的问题。

2.2 MILP的建模和求解

2.2.1 MILP第1阶段

考虑到量测量的不确定性，式（5）可改写为如下形式：

其中，y为式（5）中的辅助状态向量，其第i个分量为yi；Bi为式（5）中 B 的第 i行；分别为用以衡量量测不确定度的下界和上界；M为任意充分大的标量；bi为量测量yi对应的0／1二值变量，对于正常量测bi=0，对于不良数据bi=1。

将线性化量测不等式（6）应用于MINP，可得MILP模型如下：

MILP 模型式（7）可采用商用软件 CPLEX[30]求解，即可辨识不良数据，并可得到辅助状态向量的估计值（记为），进一步可采用BSE的第2阶段和第3阶段[27-29]得到SE的原始状态向量x的估计值，但本文不采用此法，以下给出一种更为高效的MILP第2阶段和第3阶段的建模方法。

2.2.2 MILP第2阶段

MILP的第2阶段采用如下形式的非线性变换：

利用非线性变换（8）可得到所有节点电压幅值的估计值以及所有支路两端电压相角差的估计值（一条支路可得到2个相角差估计值）。MILP第3阶段的任务是利用所有支路两端电压相角差的估计值估计除参考节点外的所有节点的相角。

2.2.3 MILP第3阶段

理论上，所有支路两端电压相角差的估计值与所有节点相角的关系为：

由于θb由MILP的第1和第2阶段的估计值计算得到，因此式（9）并不严格成立。可将由MILP的第1和第2阶段得到的θb看作含有噪声的伪量测，从而可得：

其中，τєRb为噪声矢量。

显然，式（11）可用WLS或 WLAV估计θ，这里采用WLS进行估计，可得如下估计模型：

其中，Wθ为对角权重矩阵，其对角元素取决于噪声τ的方差。

严格意义上，Wθ的对角元素应通过误差传播理论由式（8）得到，本文采用简化处理方法，令Wθ=I，则由模型（12）可得正则方程：

由于矩阵A的条件数一般较小，可直接采用正交分解法求解式（13）。模型式（13）属于二次规划QP（Quadratic Programming）问题，也可使用 CPLEX求解，即可得到所有节点相角的估计值。

综上可得MILP的计算步骤如下。

步骤1 形成矩阵B、A：形成常数雅可比矩阵B及降阶支路-节点关联矩阵A。

步骤2 求解LP：利用CPLEX求解模型式（7）。

步骤3 非线性变换：运行非线性变换式（8），得到所有节点电压幅值的估计值以及所有支路两端电压相角差的估计值（一条支路可得到2个相角差估计值）。

步骤4 求解QP：利用CPLEX求解模型式（12），得到所有节点相角的估计值。

步骤5 结束。

本文提出的MILP与MINP的理论基础相同。与MINP相比，MILP具有以下3个特点。

（1）模型仅需要求解1个MILP问题、1个非线性变换及1个QP问题，从数学上可保证获得全局最优解。

（2）算法无需非线性迭代，故不存在收敛性问题。

（3）第1阶段和第3阶段均可使用CPLEX求解，计算效率很高[30]。而与BSE相比，本文提出的方法具有良好的抗差性，可有效辨识包括杠杆点不良数据在内的多不良数据；而BSE仅可辨识单一不良数据或弱相关的多不良数据。

3 算例分析

本节以一个3节点系统及IEEE系统对所提MILP的性能进行测试，测试环境为PC机，CPU为Intel（R） Core（TM） i3 M370，主频为 2.40 GHz，内存为4.00 GB。

3.1 3节点系统算例测试

3.1.1 正确性测试

首先以图1所示的3节点系统进行测试，支路参数及量测量（均为标幺值）分别如表1和表2所示，其中节点1为参考节点。

图1 3节点系统单线图及量测配置图Fig.1 One-line diagram and measurement configuration of 3-bus system

表1 3节点系统支路参数Table 1 Branch parameters of 3-bus system

表2 3节点系统量测值Table 2 Measurements of 3-bus system

a.非线性WLS及MINP估计结果。

对图1所示的3节点系统，在如表2所示的量测集下，WLS经过4次迭代后收敛（收敛精度为10-6），其估计结果如表3所示（表中电压幅值为标幺值；为便于比较，表3同时给出了状态变量的真值）。MINP的估计结果与WLS的估计结果几乎完全一致，表中未再列出。

表3 WLS和MILP得到的状态向量估计值Table 3 Estimated state vectors by WLS and MILP

b.MILP估计结果。

采用本文提出的MILP时，第1、2阶段的估计结果如表4所示（表中，yi和分别为向量y和的第i维，分别为的估计值）。MILP第3阶段的估计结果已列于表3中。

表4 MILP第1阶段和第2阶段得到的估计结果Table 4 Estimations obtained at first and second stages of MILP

由表3可见，系统正常运行时MILP得到的状态向量估计值与WLS及MINP得到的估计值相差不大。在计算效率方面，MILP是MINP的近11倍（具体地，两者的计算耗时分别是7 ms和76 ms），MILP是WLS的4倍。而MILP的第2和第3阶段的计算效率是BSE第2和第3阶段效率的1.5倍。

3.1.2 全局寻优性测试

前文中已经指出，MILP可以从理论上确保获得全局最优解，而已有的非线性非凸形式的SE方法则无法从理论上确保获得全局最优解，本小节将给出算例证明。

仍以图1所示的3节点系统为例，考虑一组电压水平较低的情形，此时量测量如表5所示（事实上，此时系统发生了静态电压失稳现象），量测量皆为标幺值。分别使用WLS以及MILP进行估计，前者经过18次迭代后收敛（收敛精度为10-6）。状态变量真值、WLS的估计结果以及MILP的估计结果如表6所示。

表5 3节点系统另一组量测值Table 5 Another measurement set of 3-bus system

表6 WLS及MILP估计结果Table 6 Estimated state vectors by WLS and MILP

由表6可见，WLS收敛于其中一个没有实际物理意义的局部最优解；而本文提出的MILP则得到了具有实际物理意义的全局最优解，从而验证了理论分析的正确性。

3.1.3 抗差性测试

将支路1-3的电抗值缩小为原值的1／10（则节点1附近的量测量变为杠杆点）；同时在表2中的量测量Q13上增加10%的噪声以模拟不良数据，此为含杠杆点不良数据的情形。分别使用WLS+LNR（WLS带LNR辨识不良数据环节）及本文提出的MILP进行估计。

当用WLS+LNR进行估计时，首次得到的最大标准化残差对应的量测量为P2；删去P2，再次运行WLS+LNR，此时得到的最大标准化残差对应的量测量为Q31。

由此可见，对于此含有杠杆点不良数据的算例，WLS+LNR发生了残差污染和残差淹没现象，导致不良数据辨识失败。而运用本文提出的MILP则正确地辨识出不良数据Q13，从而证明了MILP可以有效抑制包括杠杆点不良数据在内的不良数据，显示了良好的抗差性。

3.2 IEEE系统算例测试

3.2.1 抗差性能测试

当噪声标准差为0.001 p.u.时，对IEEE各测试系统，由MILP得到的电压幅值和电压相角的最大估计误差（分别用‖ΔUm‖和‖Δθ‖表示）如表7所示，其中‖ΔUm‖为标幺值。

由表7可见，在所有的测试中MILP的估计结果均正确。而随着噪声标准差的减小，由MILP得到的状态向量的最大估计误差也随之减小，从而证明了MILP的估计值逼近于真值。当在以上各测试系统中随机叠加5%比例的不良数据时，MILP得到的状态向量的最大估计误差基本不变，从而证明了MILP具有良好的抗差性。

表7 MILP得到的状态变量的最大估计误差Table 7 Maximum estimation deviations of MILP

进一步在IEEE 300节点系统中，将支路1-5电抗缩小为原来的1／10（则节点1附近的量测量变为杠杆点），同时设置 4 个强相关不良数据（P1、Q1、P15和Q15），然后分别采用WLAV和MILP进行估计，估计结果如表8所示，量测量皆为标幺值。由表8可见，WLAV无法有效辨识杠杆点不良数据；而MILP则对杠杆点不良数据具有很好的抑制性，显示了良好的抗差性。

表8 IEEE 300节点系统WLAV及MILP估计结果Table 8 Estimation results of WLAV and MILP for IEEE 300-bus system

3.2.2 计算效率测试

文献[26]指出，MINP的计算效率特别低，对规模稍大系统的计算效率要求也难以满足。本小节分别在不同规模的IEEE系统上对MINP和MILP进行了多次测试，不同系统对应的量测量的个数及平均计算耗时如表9所示，MILP平均计算耗时与系统规模的关系如图2所示。由图2可见，MILP的平均计算耗时与系统规模近似成正比，故MILP适应于大规模系统的在线应用。

表9 不同IEEE测试系统中量测量的个数及平均计算耗时Table 9 Measurement quantity and average time consumption for different IEEE bus systems

图2 MILP平均计算耗时与系统规模的关系Fig.2 Relationship between average time consumption by MILP and system size

4 结语

针对已有状态估计存在的难以保证全局最优解及可能收敛困难等问题，本文提出了MILP形式的状态估计方法，该方法的优点是：（1）对包括杠杆点不良数据在内的多不良数据具有良好的抑制能力，显示了良好的抗差性；（2）从理论上可确保获得全局最优解；（3）不存在收敛性问题；（4）模型可采用商业软件CPLEX进行求解，具有很高的计算效率。