1、基于奇异值分解的同调机群识别方法 朱乔木 陈金富 段献忠 游昊 李本瑜 强电磁工程与新技术国家重点实验室(华中科技大学) 云南电力调度控制中心 摘 要： 提出一种基于奇异值分解的同调机群识别方法。该方法直接面向发电机实时功角数据, 利用奇异值分解技术提取反映机组同调性的关键信息。从能量的角度出发, 根据定义的能量贡献率指标自适应地选择奇异值数量, 构造低维权矩阵, 进而实现数据的显著降维。最终通过对权矩阵进行聚类分析实现同调机群识别。该方法具有原理简单、易于实现、计算量小等特点, 能够十分方便地实现大型复杂电网中机群的同调识别。提出了该方法的在线应用框架, 测试结果表明该方法分群速度快, 具备
2、在线应用潜力。此外, 同调识别的结果可直观地进行图形展示, 有利于电网的运行控制与分析。通过 IEEE 39 节点系统和南方电网算例验证了该方法的有效性和正确性。关键词： 电力系统; 同调机群识别; 奇异值分解; 聚类分析; 作者简介：朱乔木男, 1992 年生, 博士研究生, 研究方向为动态电力系统分析与计算, 人工智能技术在电力系统中的应用等。E-mail:作者简介：陈金富男, 1972 年生, 副教授, 研究方向为电网分析计算方法及其软件开发、电网继电保护整定计算软件开发、电力电子技术在电力系统中的应用等。E-mail:收稿日期：2016-12-01A Coherent Generato
3、rs Identification Method Based on Singular Value DecompositionZhu Qiaomu Chen Jinfu Duan Xianzhong You Hao Li Benyu State Key Laboratory of Advanced Electromagnetic Engineering and Technology Huazhong University of Science and Technology ; Yunnan Electric Power Dispatching and Controlling Center ; A
4、bstract： This paper presented a novel approach to identify coherent generators using singular value decomposition (SVD) in power systems. The key information reflecting coherency was extracted by SVD directly from the real-time power angle values provided by wide-area measurement system (WAMS) . Fro
5、m the point of energy, the weight matrix with low dimensions was constructed by referring to several main singular values, the number of which was adaptively determined by energy contribution rate defined in the paper. Coherent generators were identified by applying cluster analysis on the weight ma
6、trix. The proposed method is simple, and entails small amount of calculation, so that it is suitable for complex interconnected power grid especially. Meanwhile, this paper discussed the on-line application of the method. The results show that it has high computational speed and potential to be appl
7、ied online. Furthermore, the results of coherency identification can be expressed by graphs, which is beneficial to control and analysis of power systems. The effectiveness and correctness of the proposed method was validated by simulation results on IEEE39-bus system and China southern power grid (
8、CSG) .Keyword： Power system; coherency identification; singular value decomposition; clustering analysis; Received： 2016-12-010 引言快速、准确的同调机群识别在电力系统的动态计算、分析和控制任务中意义重大。依赖于同调机群识别的动态等值可大大降低系统规模, 减少计算量, 加快暂态稳定计算速度1。同时, 同调机群的识别结果可以揭示系统的振荡模式2,为系统的控制提供重要依据。此外, 主动解列控制也依赖于同调机群识别来确定系统中的解列点, 从而将失步机群从系统中切除, 防止电力
9、系统崩溃3-7。传统的同调机群识别方法主要基于时域仿真轨迹或线性化系统分析8,9。它们一般针对特定故障的离线仿真数据, 对系统参数依赖程度较高, 且计算代价较大10, 难以满足电力系统的实时分群要求, 在实际工程应用中受到较大限制。近年来, 广域测量系统 (Wide-Area Measurement System, WAMS) 的迅速发展和广泛应用为同调机群的识别提供了崭新的思路, 即从系统动态实时数据中挖掘机组同调信息。该技术路线不再依赖系统元件的模型和参数, 将同调机群识别从离线计算推向了在线应用, 可为系统的实时分析与控制提供依据。已有较多研究直接面向 WAMS 实测数据进行同调机群识别
10、, 例如, 文献2构建了投影寻踪模型, 通过寻找并提取反映机组同调性的关键投影方向实现同调机群识别。文献10, 11分别应用独立分量分析技术和主成分分析技术提取系统关键动态信息, 对电力系统中的机组同调性进行辨识。文献12通过连续小波变换求取小波相位差来识别同调机群。文献13通过计算各机组和母线之间的信号相关系数实现同调机群识别。但上述方法在实际应用中还存在瓶颈, 例如, 文献2, 10, 11涉及寻优过程, 计算复杂程度较高, 同调识别效率有待提高。文献12计算量较大。文献13中相关系数阈值的确定缺乏科学指导, 具有一定的主观性。本文提出一种基于奇异值分解 (Singular Value D
11、ecomposition, SVD) 的同调机群识别方法。以发电机实时功角信息为基础, 利用奇异值分解技术提取反映发电机组动态行为相似性的关键信息, 并根据定义的能量贡献率指标自适应选择奇异值的保留个数, 将维数较高的功角时间序列变换为一个维数较低的权重系数矩阵 (权矩阵) , 通过对低维的权矩阵进行聚类分析达到同调机群识别的目的。该法原理简单, 易于实现, 计算效率高, 能够方便、快速地实现复杂电网的发电机同调识别。IEEE 39 节点系统和南方电网算例仿真结果表明了该方法的有效性。1 奇异值分解的基本原理矩阵的奇异值分解是一种正交变换, 它在数据降维与压缩14,15、机械故障诊断16,17
12、、图像处理18,19、生物医学20、信号消噪21,22等诸多领域中应用广泛。基于奇异值分解的信号分析方法可以对信号进行重构, 较好地提取出特征信息, 以达到研究和分析高维数据的目的。如文献14将 SVD 用于对大批量卫星遥感数据的特征降维, 结果表明, SVD 可以获得比主成分分析更加精简和准确的特征, 减少了存储和处理量。奇异值分解的基本原理简述如下23。设矩阵 AC, rank (A) =r, 矩阵 AA 与矩阵 AA 的特征值相等, 均为 1 2 r0, r+1= r+2= n=0, 称正数 为矩阵 A 的奇异值。其中, A 为 A 的共轭转置矩阵。定理 1:设矩阵 AC, rank (
13、A) =r, 1 2 r0 是矩阵 A 的奇异值, 则存在酉矩阵 UC, VC, 分块矩阵 , 使式中, =diag ( 1, 2, , r) 。定理 2:设矩阵 AC, A 的奇异值为 1 2 r0, 令 u1, u2, , ur为相应于奇异值的左奇异向量, v 1, v2, , vr为相应于奇异值的右奇异向量, 则式 (2) 给出的形式被称为矩阵 A 的奇异值展开式, 对一个整数 kr, 略去 A的一些小的奇异值对应的项, 取矩阵 Ak, 则Ak是一个秩为 k 的 mn 阶矩阵。数学上可以严格证明, A k是在所有的秩为 k的 mn 阶矩阵中, 从 Frobenius 范数的意义上, 与矩
14、阵 A 距离最近的一个矩阵23, 即由式 (4) 和式 (5) 可知, 在一定的精度范围内, 可以用 Ak来逼近原矩阵 A, 这一原理在实际中应用非常广泛。例如, 在图像数字化技术中, 保留部分奇异值来存储图片可以达到大大压缩数据存储量, 又保持图像不失真的目的。从信息论的角度看, 信息包含于奇异值及其对应的奇异向量中, 且数值较大的奇异值对应着主要信息, 而数值较小的奇异值对应次要信息。因此, 奇异值分解可以用来挖掘高维数据中的主要信息, 这也正是本文思路的出发点。2 基于奇异值分解的同调机群识别2.1 同调机群识别问题描述电力系统的发电机同调性反映的是动态过程中系统中机组具有相似运动趋势3
15、。在关注的时间尺度 t0, 内, 如果发电机 i 和 j 的功角增量的偏差在任意时刻都不大于界定的阈值 , 则判定这两台发电机在此时间段内同调, 即式中, i (t) 、 j (t) 分别为 t 时刻发电机 i、j 的功角; i0、 j0为稳态下的功角。功角表征了发电机运动的基本特征, 包括位置、速度等信息, 实时功角可以通过 WAMS 的量测得到24。本文选取电网中多台发电机的功角作为基础数据, 进行特征降维和聚类分析, 并对其进行可视化解释, 从而确定机组的同调特性。2.2 数据预处理设获取了系统中 p 台发电机的功角量测数据 X=x1, x2, , xp, xi (i=1, 2, , p
16、) 对应第 i 台发电机, 其包含 n 个时刻上的采样值, 即 xiR, 发电机的功角时间序列为式中, n 为时间序列长度;T 为采样周期。正如式 (6) 的定义, 同调机群识别强调的是发电机在动态过程中的相似运动趋势, 因此, 在进行分群研究时, 往往更加关心发电机功角增量而非其绝对数值。故根据功角时间序列 X 计算形成功角增量时间序列 Y, 来表征各发电机功角较其稳态时刻的增量, 其计算规则为式中, Y ij为矩阵 Y 的第 i 行第 j 列元素。2.3 基于 SVD 的同调识别建模矩阵 Y 虽蕴含了机组同调信息, 但其包含了大量冗余信息。因此, 利用奇异值分解技术对矩阵 Y 进行降维,
17、提取关键信息, 可用于判别各机组的同调关系。在实际系统中, 发电机动态行为各异, 往往不存在两条变化完全一致或者严格呈比例的功角增量曲线。因此, 矩阵 Y 为列满秩矩阵, 即 rank (Y) =p。由于矩阵 Y 为实数矩阵, 依照式 (1) 可将分解结果记为式中, V 为 V 的转置矩阵。取整数 krank (Y) , 可构造近似矩阵 Yk。根据物理意义, Y k为功角增量近似时间序列, 其列向量与发电机一一对应。对Yk按列进行分块, 则任一列可表示为式中, Y k (j) 为 Yk的第 j 列。按照式 (11) , 将 ui视为基, 则 Yk的任一列 Yk (j) 都可表示为 ui的加权和
18、, 权系数为 ivji。运动特性相似的发电机对应着相近的权系数, 因此根据权系数进行归类即可实现机群的同调识别。2.4 基于 SVD 的同调识别的数学解释近似矩阵 Yk列空间的秩为 k, 基为u 1, u2, , uk。基中的向量有着特定的数学含义, 它们是一簇形状各异曲线, 通过不同的权系数线性组成功角增量曲线, 因此, 它们可以被视为功角增量曲线的基本形态 (profile) 。当形态越多时, 功角增量曲线的变化形式越多样, 变化细节越丰富。因此, k 的增加实际上就是基的扩增, 也就是 Yk列空间的拓展和升维。Y k列空间维数越高, 则对功角增量描述越精确, 反之, 则越粗糙。图 1 展
19、示了某系统中 10 台发电机功角增量曲线。由图可知, 扰动发生后, 10台发电机的功角变化规律不尽相同。当 k=1 时, Y 1的任一列都由 u1线性表示, 曲线如图 2 所示。不难发现, 各条曲线按比例同步变化, 仅在幅值上存在差异, 尽管曲线形状得到了大致刻画, 但其变化细节并未得到体现。因此, 一维空间只是刻画了全部曲线的“共性”, 而忽略了不同曲线的“个性”。为了进一步描述各曲线的变化细节, 需要更多的基向量进行组合, 因此需要将曲线放到更高维的列空间中进行表达。以图 1 中的 G1 为例, 分别取 k=1, 2, 3 时, 其功角增量曲线如图 3 所示。随着 k 增大, 近似曲线不断
20、逼近原曲线, 当 k 取 3 时, 近似曲线与原曲线差异已经非常小, 基本可以表征发电机功角的变化情况。图 1 发电机功角增量曲线 Fig.1 Curves of rotor angle increment 下载原图图 2 发电机功角增量近似曲线 Fig.2 Approximate curves of rotor angle increment 下载原图图 3 G1 功角增量近似曲线 Fig.3 Approximate curve of rotor angle increment of G1 下载原图因此, 从向量空间的角度来说, Y k列空间中的任一向量均可由 u1, u2, , uk惟一地
21、线性表示, 权系数 ivji即为 Yk (j) 在这组基下的坐标。根据坐标, 任一发电机都可映射成 k 维正交空间里的一个点。自然地, 任意两台发电机在空间中的距离则反映了两者的同调程度。2.5 k 值的确定和同调识别的步骤在低维的列空间里研究同调问题, 能给实际的分析控制带来很大的便利。当维数越低时, 原始数据矩阵的压缩程度越高, 能够减少存储空间, 节约计算资源, 简化之后的聚类分析14。本文从能量角度来自适应地确定 k 值, 用矩阵的Frobenius 范数来表示发电机功角增量时间序列的能量。由于 Frobenius 范数是酉不变范数25, 故在此基础上, 定义能量贡献率 ck以反映前
22、k 个奇异值对序列能量的贡献程度。自然地, 第 k 个奇异值对序列能量的贡献程度则可以用 ck逐差表示, 将其定义为能量增长率 ek。式中, 1kr。显然, c 1c 2c r, 也就是说, 奇异值保留越多, 则越逼近原序列的能量, 进而可以通过设定 ck的阈值来自适应地确定 k。选定 k 后, 即可根据式 (11) , 利用权系数构造权矩阵 MR, 其构造形式如式 (14) 所示。M 的维数是 kp, k 为保留奇异值的数量, 往往较小。而发电机功角时间序列X 的维数为 np, n 为采样点的数量, 往往较大。例如, 设定时间尺度为 10 s, WAMS 的刷新周期取典型值 10 ms26,
23、 k 取 5, 则 M 的规模仅为 X 的 0.5%。降维所得到的 np 阶的矩阵 M 即为同调机群识别所需要的数据。可以看出 SVD 解决了原数据维数高的问题, 实现了对原数据的特征提取。进而, 通过对权矩阵M 进行聚类分析达到同调机群识别的目的。本文采用经典的 k-means 算法12进行聚类分析。作为最经典的聚类算法之一, k-means 算法具有复杂度低、效率高、可扩展性强等优点。其基本思想是首先确定若干个初始聚类中心, 然后逐步调整这些中心, 使聚类趋于合理。值得指出的是, k-means 算法在数据含噪声、较多冗余信息、离群点或异常值时会对初始点的选择呈现出一定敏感性27。当初始点
24、选择不合理时, 算法可能出现收敛速度下降, 聚类结果容易陷入局部最优等情况28。本文中, k-means 面向的权矩阵 M 由 SVD 处理过的数据构建, SVD 本身具有的降噪、平滑和去冗余特性在很大程度上削弱了聚类环节的初值敏感性。因此, 本文采用随机化的方式选取初始点也能获得较为理想的聚类效果。此外, 虽然在 2.3 节分析推导过程中用到了 Yk, 但实际上并不需要计算 Yk, 同时也避免了求解和存储高阶方阵 U, 而只需计算奇异值及其对应的右奇异值方阵 V, 构造出权矩阵 M 即可。综上, 基于奇异值分解的同调机群识别步骤为:1) 利用 WAMS 量测数据形成发电机功角时间序列 X。2
25、) 计算得到发电机功角增量时间序列 Y。3) 对 Y 进行奇异值分解, 计算能量贡献率 ck, 根据设定的 ck阈值确定 k。4) 保留 k 个奇异值, 构造权矩阵 M。5) 对 M 进行 k-means 聚类, 得到分群结果。3 算例分析本文分别以 IEEE 39 节点系统和南方电网为例来验证本文方法的有效性。仿真工具为 PSD-BPA, 仿真时长为 10 s, 发电机功角数据采样周期为 10 ms。为了兼顾数据压缩和数据保真, 设定能量贡献率阈值为 98%。3.1 IEEE 39 节点系统算例IEEE 39 节点系统29包含 10 台发电机, 将发电机 G39 设为参考机。0.2 s 时,
26、 线路 l16-21靠近 BUS-16 处发生三相短路故障, 持续 0.2 s 后故障被切除, 故障切除后的系统保持稳定运行。系统中发电机的功角增量曲线如图 4 所示。不难看出, 当系统发生故障后, 发电机之间存在着明显的分群现象。根据功角增量曲线, 若要将发电机划分为 4 群, 则 G35 和 G36 为一群, G38 和 G39 各自单独成群, 剩余发电机为一群。图 4 算例 1 发电机功角增量曲线 Fig.4 Curves of rotor angle increment of case 1 下载原图采用本文方法进行同调识别, 对功角增量时间序列进行奇异值分解, 共得到 9个奇异值, 其
27、能量增长率和能量贡献率如图 5 所示。根据能量贡献率阈值选取前 5 个奇异值, 构造权矩阵。该矩阵维数为 510, 仅为功角增量时间序列元素量的 5, 对其进行 k-means 聚类分析, 分群结果见表 1。图 5 算例 1 奇异值能量贡献率 Fig.5 Energy contribution rate of case 1 下载原图表 1 算例 1 发电机分群结果 Tab.1 Coherency identification result of case 1 下载原表 表 1 所示的分群结果与图 4 所示的功角增量曲线是一致的, 分群结果合理。由于列空间的基是彼此正交的向量, 因此, 可以使用
28、直角坐标系对分群结果进行直观表达。在本算例中, 前三个奇异值对应的能量贡献率高达 94.56%, 已基本能够反应机组间的同调关系。本文提取前三个奇异值对应的权系数, 将其投影至三维直角坐标系中, 如图 6 所示。图中任意两点间的空间几何距离则代表了对应发电机之间的动态行为相似程度。距离越小, 相似程度越高, 反之, 则相似程度越低。例如, G35 和 G36 距离很近, 而与其他机组空间距离较大, 说明这两台发电机同调程度较高, 而且其运动行为明显区别于其他机组。这样的分析结果, 与功角增量曲线一致。此外, 根据能量贡献率不同, 也可以将权系数投影至二维平面或者一维区间进行分析。图 6 算例
29、1 同调机群三维投影 Fig.6 3-D projection of coherency identification of case 1 下载原图3.2 南方电网算例采用 2016 年南方电网夏季大方式作为算例, 验证本文所提方法在大系统中的有效性。2016 年, 云南电网与主网 (广东、广西、贵州、海南) 仅通过直流互联, 二者实现异步运行。对于南方电网实际运行来说, 同调机群识别将为振荡模式分析、稳控措施的制定等提供重要依据。由于该算例规模较大, 机组众多, 本文篇幅有限, 故挑选 70 台发电机组作为研究对象, 见表 2。表 2 研究对象 Tab.2 Study objects 下载原
30、表 故障设置为广东区域内 DH 线单相短路, 500 k V HD 站单相中开关拒动。故障持续 15 周波, 同串线路跳开后该故障才得到切除。同时, 该故障将导致 8 回直流短时换相失败, 给系统带来较大扰动。作为两个同步电网, 主网和云南区域内机组功角增量曲线分别如图 7a、图 7b 所示 (参考机分别取 SJC、LH) 。由图可知, 由于多回直流换相失败, 云南电网与主网之间的送电通道被短时切断, 两区域内的功率平衡被打破, 各发电机组功角出现较为剧烈的变化。图 7 算例 2 机组功角增量曲线 Fig.7 Curves of rotor angle increment of case 2
31、下载原图对云南电网和主网机组的功角增量时间序列进行奇异值分解, 分别得到 16 和54 个奇异值, 其 ek和 ck见表 3, 由于数值差别较大, 为了清楚表示, 仅列出前5 个奇异值。根据设定的能量贡献率阈值, 主网和云南分别选取前 2 个和前 3个奇异值构造权矩阵, 权矩阵的元素个数仅分别为功角增量时间序列的 2和3。表 3 能量贡献率 Tab.3 Energy contribution ratio 下载原表 由表 3 中可知, 对于主网机组来说, 前两个奇异值的能量贡献率已逾 98%, 可认为前两个奇异值对应的权系数即可以基本反映系统同调模式。提取前两个奇异值对应的权系数, 将 54 台发电机投影至二维平面, 结果如图 8 所示, 图中标识出了将机组分作两群时的同调识别结果。由图可知, 单利用 x 坐标即可基本揭示该故障激发的贵 (州) 广 (东) 振荡模式, 即贵州机组与广东机组相对振
