简介这份资源面向电力系统分析与智能电网方向的学习者和研究人员聚焦谱聚类算法在电力网络分区中的实现帮助解决复杂电网结构下节点划分与区域划分的建模问题。压缩包共2个文件均为MATLAB源码.m整体约1KB体量轻便便于直接阅读与二次修改。内容围绕IEEE 30节点标准系统展开涉及电气距离计算、相似度与权重矩阵构建、拉普拉斯矩阵特征分解以及特征向量聚类等关键环节可帮助读者理解从电气距离到分区结果的完整算法链路。已有503人学习下载说明该案例在电力系统分区研究中具有一定参考价值。对于需要复现谱聚类分区流程、验证电气距离指标效果或开展课程设计、论文实验的读者这份源码可作为算法实现与调试的起点也便于在此基础上调整参数、对比不同分区方案。1. 电力系统分区遇上谱聚类从电气距离到集群划分的完整落地做电网分析的人大多有过这样的经历拿到一张几十甚至上百节点的网络拓扑想按电气特性把它切成几个相对独立的区域用传统聚类方法跑出来的结果却总是不尽如人意——要么某个分区里塞进了电气上根本不亲密的节点要么分区边界来回震荡换个初始值结果就变了个样。电力系统分区这件事核心诉求其实很朴素让同一分区内的节点电气联系紧密不同分区之间联系稀疏这样后续做无功优化、状态估计、电压控制才能各管各的、互不干扰。谱聚类算法之所以在这个场景里被反复提起是因为它把分区问题转化成了图割问题用电气距离构造权重矩阵再借助拉普拉斯矩阵的特征分解完成降维和切分天然适合处理非凸、形状不规则的簇结构。这份资源围绕谱聚类分区和电气距离两个关键词展开适合正在做电网分区、无功电压控制分区或者对谱聚类落地感兴趣的从业者新手能照着流程走通熟手能从中看到参数边界和工程取舍。2. 谱聚类做电力系统分区原理拆解与电气距离矩阵构造2.1 为什么电力系统分区偏偏选中了谱聚类传统分区方法大致分几类基于灵敏度矩阵的逐次分裂、基于模块度的社区发现、以及各种 K-means 变体。它们各有各的问题。灵敏度法对初始分区敏感分裂顺序不同结果差异明显模块度方法在电网这种带权图上容易陷入分辨率极限小分区被合并、大分区被拆散K-means 则默认簇是凸的、各向同性的而电力系统的电气关联结构往往是非凸的强行用欧氏距离聚类等于假设两个节点电气距离近就一定要分到一起忽略了网络拓扑的连通性约束。谱聚类的思路不一样。它先把网络建成一张带权无向图 G(V, E)V 是节点母线或区域E 是边边的权重 w_ij 反映节点 i 和 j 之间的电气关联强度。然后构造拉普拉斯矩阵对矩阵做特征分解取前 k 个特征向量组成新的特征空间最后在这个低维空间里做一次 K-means 或者离散化切分。这样做的好处是特征分解过程自动捕捉了图的全局连通结构分区结果对初始值不敏感而且能处理任意形状的簇。对于电力系统分区来说这意味着分区边界更符合电气上的自然分界不会出现“飞地”式的孤立节点被硬塞进某个分区的情况。选型上还有一个现实考量谱聚类对相似度矩阵的构造方式很敏感而电力系统恰好有明确的电气距离定义可以构造出物理意义清晰的权重矩阵。这一点比纯数据驱动的聚类方法有优势——你不是在盲目调参而是在用电气距离这个有物理含义的量来引导分区。2.2 电气距离怎么算从阻抗矩阵到权重矩阵电气距离是这份资源的核心概念之一。在电力系统里两个节点之间的电气距离不是地理距离也不是简单的线路阻抗而是反映它们电压耦合程度的量。常见做法是基于节点阻抗矩阵 Z_bus 来构造。假设系统有 n 个节点节点阻抗矩阵 Z 是 n×n 的方阵其中 Z_ij 表示在节点 j 注入单位电流时节点 i 的电压响应。电气距离的一种常用定义是d_ij Z_ii Z_jj - 2·Z_ij这个量在电气上等价于节点 i 和 j 之间的等效电气距离数值越小表示两个节点电气耦合越紧。注意 Z_ij 通常取阻抗矩阵的模或者实部具体取法要看应用场景——做无功电压分区一般取模值做有功相关分析可能取电抗部分。有了电气距离矩阵 D下一步是把它转成谱聚类需要的相似度矩阵 W。因为谱聚类要求权重越大表示越相似而电气距离越大表示越远所以需要一个转换。常见做法是高斯核w_ij exp(-d_ij² / (2σ²))其中 σ 是尺度参数控制相似度衰减速度。σ 取太小权重矩阵会退化成近似对角阵分区结果碎片化σ 取太大所有权重趋近于 1分区失去区分度。工程上一般用电气距离矩阵的中位数或者均值来标定 σ也可以做几次试算看分区数的稳定性。下面是一段构造电气距离矩阵和相似度矩阵的 Python 代码假设你已经有了节点阻抗矩阵 Zimport numpy as np from scipy.linalg import eigh def build_electrical_distance(Z): 从节点阻抗矩阵构造电气距离矩阵 Z: n x n 复数矩阵节点阻抗矩阵 返回: n x n 实数电气距离矩阵 n Z.shape[0] # 取阻抗矩阵的模作为电气距离计算基础 Z_abs np.abs(Z) D np.zeros((n, n)) for i in range(n): for j in range(n): # 电气距离定义: d_ij Z_ii Z_jj - 2*Z_ij D[i, j] Z_abs[i, i] Z_abs[j, j] - 2 * Z_abs[i, j] # 对称化处理消除数值误差 D (D D.T) / 2 np.fill_diagonal(D, 0) return D def build_similarity(D, sigmaNone): 用高斯核把电气距离转成相似度矩阵 D: 电气距离矩阵 sigma: 尺度参数默认取距离矩阵的中位数 if sigma is None: # 取上三角非零元素的中位数作为尺度 triu_vals D[np.triu_indices_from(D, k1)] sigma np.median(triu_vals) W np.exp(-D**2 / (2 * sigma**2)) np.fill_diagonal(W, 0) # 自相似置零避免自环影响 return W这段代码里有两个关键参数需要留意。一个是电气距离的定义方式上面用的是 Z_ii Z_jj - 2·Z_ij 这种形式实际工程中也有用 1/|Z_ij| 或者 |Z_ii - Z_ij| 的变体选哪种取决于你的分区目标——如果关注的是节点间的电压相互影响前者更合适如果关注的是节点对某参考点的耦合强度后者可能更直接。另一个是 σ 的取值代码里默认用中位数这是一个比较稳健的起点但如果你的系统节点数很少比如少于 20 个中位数可能偏大需要手动调小。提示节点阻抗矩阵 Z_bus 的求取本身就是一个独立步骤通常用潮流计算或者短路计算程序生成。如果你手头只有导纳矩阵 Y_bus可以先求逆得到 Z_bus但要注意 Y_bus 可能是奇异的比如存在孤立节点求逆前需要处理。2.3 拉普拉斯矩阵与特征分解分区是怎么“切”出来的有了相似度矩阵 W接下来构造拉普拉斯矩阵。谱聚类里常用的有两种非归一化拉普拉斯 L D_deg - W以及归一化拉普拉斯 L_sym I - D_deg^(-1/2) · W · D_deg^(-1/2)。电力系统分区一般用归一化版本因为各节点的度连接强度差异可能很大归一化能避免高度节点主导分区结果。其中 D_deg 是度矩阵对角元 D_deg[i,i] Σ_j w_ij。然后对 L_sym 做特征分解取前 k 个最小特征值对应的特征向量排除接近零的第一个特征向量它对应的是平凡解组成 n×k 的特征矩阵 U。把 U 的每一行看作一个样本做一次 K-means 聚类得到的标签就是每个节点的分区归属。k 的选取是另一个需要斟酌的点。理论上可以用特征值间隙eigengap来估计把特征值从小到大排列找到相邻特征值差值最大的位置那个位置对应的序号就是比较合适的 k。但在电力系统里分区数往往还受运行约束影响——比如一个省级电网可能要求分成 3 到 5 个无功电压控制区这时候 k 就不是纯数据驱动的了需要结合工程需求来定。def spectral_partition(W, k): 对相似度矩阵做谱聚类分区 W: n x n 相似度矩阵 k: 分区数 返回: 每个节点的分区标签 n W.shape[0] # 度矩阵 deg np.sum(W, axis1) deg_inv_sqrt np.where(deg 0, 1.0 / np.sqrt(deg), 0) D_deg_inv_sqrt np.diag(deg_inv_sqrt) # 归一化拉普拉斯矩阵 L_sym I - D^{-1/2} W D^{-1/2} L_sym np.eye(n) - D_deg_inv_sqrt W D_deg_inv_sqrt # 特征分解取前 k 个最小特征值对应的特征向量 eigvals, eigvecs eigh(L_sym) # 跳过第一个接近零的特征向量 U eigvecs[:, 1:k1] # 对特征向量矩阵的行做归一化 norms np.linalg.norm(U, axis1, keepdimsTrue) norms[norms 0] 1 U_norm U / norms # K-means 聚类 from sklearn.cluster import KMeans kmeans KMeans(n_clustersk, n_init10, random_state42) labels kmeans.fit_predict(U_norm) return labels这段代码里eigh 用于对称矩阵的特征分解比通用的 eig 更稳定也更快。取特征向量时跳过了第一个因为归一化拉普拉斯的最小特征值总是 0对应的特征向量是常向量不携带分区信息。U 的行归一化是为了让 K-means 在单位球面上工作这是归一化谱聚类的标准做法能提升数值稳定性。K-means 的 n_init 设成 10 是常见做法多跑几次取最优避免陷入局部最优。random_state 固定是为了结果可复现实际工程中如果对分区稳定性有要求可以多跑几组随机种子看标签一致性。3. 从阻抗矩阵到分区标签一套可复现的实操流程3.1 数据准备你需要哪些输入文件跑通这套分区流程最少需要两样东西网络拓扑参数和节点阻抗矩阵。拓扑参数包括支路连接关系、支路阻抗、变压器变比等这些通常来自 PSS/E、BPA 或者 MATPOWER 的算例文件。如果你用的是 MATPOWER 格式的 case 文件可以用 Python 的 pandapower 或者 pypower 库直接读取。节点阻抗矩阵 Z_bus 有两种获取途径。一是从潮流程序里导出很多商用软件都支持输出 Z_bus 或者短路计算用的阻抗矩阵。二是自己从 Y_bus 求逆前提是 Y_bus 非奇异。对于含变压器的系统注意变压器模型要正确反映到 Y_bus 里否则求出来的 Z_bus 电气距离会失真。下面是一个从 MATPOWER case 文件读取数据并构造 Y_bus 的示例import numpy as np from scipy.io import loadmat def load_case_data(case_file): 读取 MATPOWER 格式的 case 文件 返回: 支路数据、节点数据 data loadmat(case_file) # MATPOWER 的 branch 矩阵列定义: # 0: from_bus, 1: to_bus, 2: r, 3: x, 4: b, ... branch data[branch] # bus 矩阵列定义: # 0: bus_id, 1: type, 2: Pd, 3: Qd, ... bus data[bus] return branch, bus def build_ybus(branch, bus, baseMVA100): 从支路数据构造节点导纳矩阵 branch: 支路矩阵 bus: 节点矩阵 baseMVA: 基准容量 n bus.shape[0] Y np.zeros((n, n), dtypecomplex) # 建立 bus_id 到索引的映射 bus_ids bus[:, 0].astype(int) id_to_idx {bid: i for i, bid in enumerate(bus_ids)} for row in branch: f id_to_idx[int(row[0])] t id_to_idx[int(row[1])] r, x, b row[2], row[3], row[4] z r 1j * x y 1.0 / z if z ! 0 else 0 # 串联支路导纳 Y[f, f] y Y[t, t] y Y[f, t] - y Y[t, f] - y # 并联充电电容 Y[f, f] 1j * b / 2 Y[t, t] 1j * b / 2 return Y这段代码里build_ybus 只处理了最基本的串联支路和并联电容没有考虑变压器变比。如果你的系统里有变压器需要在支路循环里额外处理变比——通常做法是把变比折算到导纳矩阵的对应元素上。另外baseMVA 参数在标幺值系统里用于容量基准转换如果你的支路参数已经是标幺值这个参数不影响 Y_bus 的构造但后续如果要做功率相关的计算会用到。拿到 Y_bus 之后求逆得到 Z_bus# 检查 Y_bus 是否可逆 if np.linalg.matrix_rank(Y) Y.shape[0]: print(Y_bus 奇异可能存在孤立节点需要先处理) else: Z np.linalg.inv(Y) D build_electrical_distance(Z) W build_similarity(D) labels spectral_partition(W, k4) print(分区标签:, labels)这里有一个容易翻车的点Y_bus 奇异。电力系统里如果存在没有支路连接的孤立节点或者某个节点的对地导纳为零Y_bus 就会出现零特征值求逆直接报错。处理办法是先识别孤立节点把它们从系统里剔除或者给这些节点加一个很小的对地导纳比如 1e-6让矩阵满秩。后者会引入微小误差但工程上通常可以接受。3.2 参数调优σ 和 k 怎么定σ 和 k 是这套流程里最需要动手调的两个参数。σ 控制相似度矩阵的稀疏程度k 决定分区数。两者不是独立的——σ 变了特征值分布会变最优 k 也可能跟着变。我一般会按这个顺序来调先固定一个初始 k比如根据系统规模估一个50 节点以下取 3 到 4100 节点以上取 5 到 8然后扫一遍 σ 的取值看分区结果的稳定性。具体做法是取 σ 为电气距离矩阵的 0.5 倍中位数、1 倍中位数、2 倍中位数、5 倍中位数分别跑分区看标签的调整兰德指数ARI或者归一化互信息NMI。如果几个 σ 下结果高度一致说明分区结构稳健如果差异很大说明 σ 需要更精细地调。k 的选取可以结合特征值间隙def estimate_k(W, max_k10): 用特征值间隙估计合适的分区数 n W.shape[0] deg np.sum(W, axis1) deg_inv_sqrt np.where(deg 0, 1.0 / np.sqrt(deg), 0) D_deg_inv_sqrt np.diag(deg_inv_sqrt) L_sym np.eye(n) - D_deg_inv_sqrt W D_deg_inv_sqrt eigvals np.sort(np.linalg.eigvalsh(L_sym)) # 计算相邻特征值差值 gaps np.diff(eigvals[:max_k1]) # 最大间隙对应的位置 k_est np.argmax(gaps) 1 return k_est, eigvals[:max_k1]特征值间隙法的逻辑是如果图里有 k 个自然分簇那么拉普拉斯矩阵的前 k 个特征值会明显小于后面的第 k 个和第 k1 个之间会出现一个较大的间隙。但实际电网里这个间隙往往不那么明显所以 k_est 只能作为参考最终还是要结合运行需求来定。注意σ 取太小的时候相似度矩阵会变得非常稀疏拉普拉斯矩阵可能出现多个接近零的特征值导致特征值间隙法失效。这时候分区结果会碎片化出现很多只有一两个节点的分区。遇到这种情况先把 σ 调大再重新估计 k。3.3 结果验证分区质量怎么评估分区跑出来之后不能只看标签就完事得验证分区质量。常用的指标有几个一是分区内电气距离均值。对每个分区计算分区内所有节点对的电气距离平均值这个值越小说明分区内节点电气耦合越紧。二是分区间电气距离均值越大越好。三是分区规模均衡度如果某个分区节点数特别少或者特别多可能需要调整 k 或者 σ。下面这段代码计算分区内和分区间的电气距离统计def evaluate_partition(D, labels): 评估分区质量 D: 电气距离矩阵 labels: 分区标签 unique_labels np.unique(labels) intra_dists [] inter_dists [] for li in unique_labels: idx_i np.where(labels li)[0] # 分区内距离 if len(idx_i) 1: intra D[np.ix_(idx_i, idx_i)] intra_vals intra[np.triu_indices_from(intra, k1)] intra_dists.append(np.mean(intra_vals)) # 分区间距离 for lj in unique_labels: if lj li: continue idx_j np.where(labels lj)[0] inter D[np.ix_(idx_i, idx_j)] inter_dists.append(np.mean(inter)) print(f分区内平均电气距离: {np.mean(intra_dists):.4f}) print(f分区间平均电气距离: {np.mean(inter_dists):.4f}) print(f分区规模: {[np.sum(labels l) for l in unique_labels]}) return np.mean(intra_dists), np.mean(inter_dists)这个评估函数会输出三个关键数字分区内平均电气距离、分区间平均电气距离、以及每个分区的节点数。理想情况下分区内距离应该明显小于分区间距离比值越大说明分区效果越好。如果两者接近说明分区没有把电气上紧密的节点有效聚在一起需要回头检查 σ 和 k 的设置。分区规模也需要看一眼。如果出现某个分区只有一两个节点而另一个分区占了系统的一大半这种分区在实际运行中很难用——要么是 σ 太小导致过度切分要么是 k 设得太大。我一般会要求每个分区的节点数不少于总节点数的 10%低于这个比例就认为分区不合理。4. 谱聚类分区避坑指南五个真实翻车场景4.1 现象分区结果每次跑都不一样原因K-means 的随机初始化没有固定或者 σ 取值处于临界区域相似度矩阵的微小数值扰动会导致特征向量方向翻转。解决固定 random_state并且在 K-means 里把 n_init 设大一点比如 20 或 50。如果还是不稳定说明 σ 需要调整——把 σ 往中位数的 1.5 倍或 2 倍方向调让相似度矩阵更平滑特征值分布更稳定。另外特征分解用 eigh 而不是 eig前者对对称矩阵的数值处理更稳健。4.2 现象某个分区里出现了地理上不相邻、电气上也看不出关联的节点原因电气距离矩阵构造时用了错误的 Z_bus。常见情况是 Y_bus 里变压器的变比没有正确折算导致某些节点的自阻抗和互阻抗计算错误电气距离失真。解决检查变压器支路的处理逻辑。变压器在 Y_bus 里的贡献不是简单的串联导纳需要根据变比在对应节点上做修正。如果用的是 MATPOWER 的 case 文件注意 branch 矩阵里的 ratio 列变比不为零的支路要特殊处理。验证方法是拿几个已知电气耦合紧密的节点对看它们的电气距离是否合理。4.3 现象特征值间隙法给出的 k 和实际需求差很远原因电网的拉普拉斯矩阵特征值分布往往没有明显的间隙尤其是节点数多、连接密集的系统。特征值间隙法在这种场景下不可靠。解决不要迷信特征值间隙。把 k 的选取和工程需求绑定——比如无功电压分区通常按变电站数量或者控制区域来定直接指定 k 的范围在这个范围内扫一遍用分区内/分区间电气距离比值来选最优 k。特征值间隙只作为参考不作为唯一依据。4.4 现象分区内电气距离和分区间电气距离差不多分区效果不明显原因σ 取得太大相似度矩阵趋近于全 1 矩阵拉普拉斯矩阵的特征向量携带的分区信息被淹没。或者电气距离矩阵本身区分度不够比如系统里所有节点的电气距离都差不多。解决先把 σ 调小从中位数的 0.3 倍开始试。如果调小 σ 后分区内距离明显下降说明之前 σ 偏大。如果调小 σ 后分区碎片化说明系统本身的电气距离分布就比较均匀这时候可能需要换一种相似度构造方式比如用 k 近邻图代替全连接高斯核只保留每个节点的前 k 个最近邻的边。4.5 现象Y_bus 求逆报错或者结果全是 NaN原因Y_bus 奇异通常是因为系统里有孤立节点或者某个节点的对地导纳为零导致矩阵秩亏。另外如果支路参数里有零阻抗支路理想变压器或者母线联络开关Y_bus 也会出现无穷大元素。解决先检查 Y_bus 的秩找出秩亏的维度。对于孤立节点直接从系统里剔除分区完成后再映射回去。对于零阻抗支路用一个小阻抗代替比如 1e-4 标幺值避免数值溢出。求逆之前先做条件数检查条件数大于 1e12 的话求逆结果不可信需要改用伪逆或者正则化求逆。5. 进阶技巧把谱聚类分区嵌入无功电压控制流程分区本身不是目的目的是为后续控制服务。在无功电压控制里分区结果直接决定了哪些节点归哪个控制器管。一个常见的进阶用法是把谱聚类分区和灵敏度分析结合起来先用谱聚类得到初始分区再在每个分区内计算无功电压灵敏度矩阵根据灵敏度符号和大小微调分区边界。具体做法是对每个分区计算分区内所有节点的电压-无功灵敏度 ∂V/∂Q。如果某个节点对分区内其他节点的灵敏度均值低于阈值说明它和这个分区的电气耦合不够紧考虑把它移到相邻分区。移动之后重新计算分区质量指标如果分区内电气距离下降、分区间距离上升就保留这次调整否则回退。这个微调过程可以用下面这段代码实现def refine_partition(D, labels, Z, threshold0.1): 基于灵敏度微调分区边界 D: 电气距离矩阵 labels: 初始分区标签 Z: 节点阻抗矩阵 threshold: 灵敏度阈值 n len(labels) labels_refined labels.copy() unique_labels np.unique(labels) for i in range(n): current_label labels_refined[i] # 计算节点 i 对当前分区内其他节点的平均灵敏度 same_partition np.where(labels_refined current_label)[0] same_partition same_partition[same_partition ! i] if len(same_partition) 0: continue # 用阻抗矩阵的模作为灵敏度代理 sens_current np.mean(np.abs(Z[i, same_partition])) # 找灵敏度更高的相邻分区 best_label current_label best_sens sens_current for lj in unique_labels: if lj current_label: continue other_partition np.where(labels_refined lj)[0] sens_other np.mean(np.abs(Z[i, other_partition])) if sens_other best_sens * (1 threshold): best_sens sens_other best_label lj if best_label ! current_label: labels_refined[i] best_label return labels_refined这段代码的逻辑是对每个节点比较它对当前分区内节点的平均灵敏度和对相邻分区内节点的平均灵敏度。如果某个相邻分区的灵敏度明显更高超过 threshold 设定的比例就把这个节点移过去。threshold 设成 0.1 表示需要高出 10% 才触发移动设太小会导致频繁调整设太大则微调效果不明显。微调之后用第 3 章里的 evaluate_partition 重新评估分区质量。如果分区内平均电气距离下降了说明微调有效。如果反而上升说明 threshold 设得不合适需要调大。还有一个实用技巧是分区编号的稳定化。谱聚类跑出来的标签是任意编号的这次跑可能分区 1 是左上角那片下次跑可能变成右下角。如果要把分区结果接入实际控制系统编号不稳定会带来麻烦。解决办法是每次分区完成后按分区内节点编号的最小值或者分区质心的电气距离排序来重新编号保证同一片区域在不同次运行中拿到相同的标签。我自己的习惯是每次跑完分区先把标签按分区内节点数从大到小排序节点数相同的按分区内最小节点编号排序。这样即使 σ 或 k 有微小变化只要分区结构大体不变编号就是稳定的。这个习惯帮我省了很多次在控制程序里追分区编号的麻烦。提示分区微调之后一定要重新验证分区规模均衡度。灵敏度微调有可能把大量节点吸到某个分区里导致分区规模失衡。如果出现这种情况给每个分区的节点数设一个上限超过上限的移动请求直接拒绝。希望这套从电气距离构造到谱聚类分区、再到灵敏度微调的流程能帮你在电力系统分区这件事上少走几个弯路。本文还有配套的精品资源点击获取