在配电网规划这个圈子里“旋转中心线距离加权交替定位算法”乍一听像某个数学竞赛题但它其实是中压配电网供电单元划分里相当实用的一招。我最近刚刚完整复现了这篇《考虑负荷特性互补及供电单元划分的中压配电网实用……》论文里的核心算法下面把整个程序复现的思考过程、代码实现和踩坑记录都整理出来。如果你也在研究配电网负荷分区、馈线供电范围优化或者想借鉴一种比K-means更适合带状区域的聚类模型这篇文章应该能帮你省掉不少摸索时间。需要先说清楚的是论文完整的模型里除了空间划分还包含网络接线、馈线容量这些偏工程化的约束。我复现的重点放在两个核心点上一是负荷点在空间上的供电单元划分二是负荷特性互补的量化处理。整个算法的名字已经透露出它的思路——“旋转中心线”表示每个供电单元用一条直线段代表其走向“距离加权”表示划分时每个负荷点按自身负荷大小影响决策“交替定位”则表示通过固定一方优化另一方的方式迭代求解。理解了这三个关键词后面的代码就不难读懂了。1. 这个算法到底在解决什么问题1.1 供电单元划分是配电网规划的底层工作中压配电网规划里有一个绕不开的环节把一片区域内大量分散的负荷点划分成若干个供电单元每个单元将来由某条馈线或某个变电站的供电范围来覆盖。这个工作如果靠人工在地图上画常常费时费力画出来的方案也很难说清哪里最优更别说应对负荷增长后的调整。供电单元划分的难点在于它既要让每个单元在地理上尽量紧凑连续又要让各单元负荷总量比较均衡避免有的馈线重载、有的馈线闲置。如果我们把负荷点当作一组带坐标和负荷数值的点那么这本质上就是一个空间聚类问题。K-means之类的方法确实可以划分但K-means默认每个簇的中心是一个点而配电网的馈线走向通常是带状延伸的用一个质心点代表一个供电单元无法体现“沿线路走廊发展”的特征。这也是为什么论文里采用“中心线”而不是“中心点”来代表供电单元。一条直线段天然能表达馈线的主干走向负荷点距离某条中心线越近说明接入对应馈线越方便网损也往往更小。把聚类中心从点换成线看起来只是几何形状改变背后的优化目标和迭代方式都要重新设计。1.2 负荷特性互补为什么比单纯看距离更重要空间距离近只是划分的一方面。实际当中同一个供电单元里的负荷如果全是同一种类型比如全是工业园区负荷白天用电很大、夜间负荷很小那么给这个单元供电的线路和变压器都要按白天的尖峰配置设备利用率其实很低。反过来如果一个供电单元里既有白天高峰的工业负荷又有晚上高峰的居民负荷两者联合在一起整体负荷曲线会平滑很多需要配置的变电容量也可以按较低的最大负荷来设计。这就是“负荷特性互补”的价值。论文把每个负荷点看作带有一条日负荷曲线的对象曲线形状不同彼此之间就能形成错峰互补。在划分供电单元时如果能把曲线相关性较高、峰谷时段接近的负荷尽量分到不同单元或者反过来说让每个单元内部呈现多种时间特性组合那么整体负荷曲线就会更平稳设备利用率能明显提升。算法里处理这一目标的方式是在空间距离代价的基础上增加一个特性互补惩罚项。每个负荷点不仅仅用坐标和大小参与划分还带一条24小时或若干时段的归一化负荷曲线当某个点准备划入某个单元时程序会计算“加入后该单元综合曲线波动性增加了多少”。增量越小说明这个点的曲线与单元内其他点越互补这个点就越应该进入这个单元。1.3 从论文提炼出的核心优化模型我复现时把问题抽象成这样一个模型有N个负荷点每个点的属性包括平面坐标x_i, y_i、负荷量P_i、以及负荷特性曲线向量C_i。目标是把N个点分成K个供电单元每个单元对应一条中心线中心线用角度θ和偏移ρ来表示直线方程是x乘cosθ加y乘sinθ等于ρ。目标函数由两部分构成。第一部分是空间代价每个点到所属单元中心线的垂直距离的平方乘以该点负荷量再加总第二部分是特性互补代价每个单元内部综合负荷曲线的方差加总。两部分通过一个权重系数λ来平衡。空间代价越小划分越符合地理和网损要求特性代价越小单元内负荷越互补。这里要说明一下论文里可能还有更严谨的数学表达和约束条件但算法求解框架基本就是“分配—更新—再分配—再更新”的交替迭代。下面一节详细拆解这个交替过程。2. 旋转中心线距离加权交替定位算法原理解读2.1 把供电单元划分看作“线聚类”问题接触过K-means的读者可以很快理解这个算法的骨架。K-means是给定K个点中心每个样本归属到距离最近的中心然后重新计算中心迭代直到稳定。这里的“旋转中心线距离加权交替定位算法”做的是同一件事只不过把“点中心”换成了“线中心”。为什么要用线前面已经提到供电单元不是各个方向均匀扩张的圆形区域而是沿线路走向延伸的狭长区域。如果用电网拓扑的语言描述一个供电单元的主干线路基本可以用一条折线近似简化成直线段后负荷点到主干线的垂直距离就代表了接入线路的代价。相比点到中心的欧氏距离点到线的垂直距离更能刻画“沿线供电”的成本。不过把中心从点换成线会带来一个新问题点的中心可以用均值直接求线的中心角度和偏移却不能用简单的平均值表达。论文里提出的“旋转中心线”方法给出了一个解析答案每个单元内所有负荷点按负荷量加权后最优中心线是过加权质心、方向与数据主方向一致的直线。主方向可以用加权协方差矩阵的特征值分解得到这一步在程序里其实只需要几行numpy代码。2.2 “旋转中心线”的几何含义与数学表达先假设我们已经把负荷点暂时分配给某个单元需要根据单元内所有点求一条最佳中心线。最小化目标是这些点到直线的负荷加权垂直距离平方和。给定一条直线法向量为n偏移为ρ。点到直线的垂直距离是x向量点乘法向量再减去ρ的绝对值。我们的目标是选择法向量n和偏移ρ使得各点到该直线的加权距离平方和最小。这个问题可以分解成两步最优偏移ρ等于所有点按负荷加权的质心在法向量方向上的投影也就是说直线必须经过加权质心最优法向量n指向数据方差最小的方向。为什么是最小方差方向因为点到直线的距离平方可以看成把点投影到法向量方向上得到的坐标与ρ的差。为了让整体投影最集中需要找到方差最小的投影方向。加权协方差矩阵的特征值分解中最小特征值对应的特征向量就是我们要的法向量。实际程序里也可以先求最大特征值对应的“主方向”然后取与之垂直的方向作为法向量这等价于先确定馈线走向再求中心线。我更喜欢“方向向量”的写法最大特征值对应的特征向量v就是馈线的主走向法向量n取v旋转90度偏移ρ等于质心点乘n。这样中心线的角度就是主方向的角度再调整90度程序里用arctan2可以方便地算出。2.3 “距离加权”与“交替定位”的迭代逻辑接下来是完整迭代框架。第一步给每条中心线随机或启发式初始化一个角度和偏移比如可以从全局数据的PCA主方向附近随机生成K组参数。第二步是分配对所有负荷点计算它到K条中心线的加权距离加权距离等于垂直距离乘以负荷量平方根或负荷量本身具体看目标函数定义我把负荷量直接乘在距离平方上分配时每个点选择加权距离最小的中心线。第三步是更新每个负荷点有了新的归属后对每个单元重新计算加权协方差矩阵求出新的中心线参数。第四步是判断收敛计算当前目标函数值与上一轮比较如果变化小于阈值就停止否则回到第二步。这个“固定中心线分配负荷点固定分配更新中心线”的结构就是“交替定位”的含义。它很像坐标下降法或EM算法里交替估计隐藏变量和参数的思路。每一次分配步骤保证目标函数不增加每一次更新中心线也保证目标函数不增加所以整个迭代过程是单调下降的至少不会劣化。值得注意的是这里的目标函数只包含空间距离部分时交替迭代是有保证收敛性质的。一旦加入负荷特性互补代价分配步骤中每个点分配到一个单元会直接影响单元综合曲线的方差而更新中心线时不涉及特性项所以整体仍然是在交替优化两个变量块不过目标函数不再是纯凸问题收敛到全局最优没有理论保证只能通过多次随机初始化来逼近。2.4 目标函数与收敛条件怎么设定在程序里我定义目标函数为J等于空间代价加λ乘特性代价。空间代价很好算每个点到中心线垂直距离的平方乘以该点负荷量再对所有点求和。特性代价稍微复杂先把每个单元内所有点的负荷曲线按负荷量加权平均得到单元综合曲线然后计算这个综合曲线在一天内的方差把所有单元的方差求和。曲线方差越小峰谷波动越平滑。收敛条件我一般看两处一是目标函数的变化量小于1e-6二是所有负荷点的归属标签与上一轮完全一致。实际跑下来只要数据规模不大通常几十轮就能稳定。如果λ设置得比较大特性代价会剧烈争夺分配权目标函数可能出现小幅度振荡这时候需要把收敛条件放松一些或者对λ做衰减。这里有个容易被忽略的代码细节中心线参数里角度θ和偏移ρ并不是一一对应的同一根直线可以用θ加π再加负ρ表示。如果不做归一化处理迭代过程中可能因为参数表达的跳变导致目标函数明明没变代码却误判为不收敛。我实现的统一做法是把θ归一化到0到π区间同时调整ρ的正负号确保每个单元的中心线只有一种表示方式。3. 程序复现用Python实现整个流程3.1 环境准备与数据结构设计复现这个算法不需要复杂框架我用的环境是Python 3.10加numpy和matplotlib核心逻辑全在numpy里可视化只用matplotlib画划分效果。如果只是复制算法不画图的话numpy一个包就够了。数据结构上我用一个N乘2二维数组保存所有负荷点坐标用一个长度N的数组保存负荷量。负荷特性曲线不能简单用一维数组因为每条曲线是24小时数据所以用N乘24二维数组保存。为了让不同量纲的数据能放到一起比较我会对坐标做标准化并把每条负荷曲线的峰值归一到1。这样距离代价和特性代价的量级通常都落在1到100之间λ调整起来比较顺手。“点、线、曲线”三个对象的关系在代码里其实不需要定义类直接用数组和字典就能表达。我习惯把状态集中在一个字典里包含坐标数组、负荷量数组、曲线数组、当前中心线参数数组以及每个点的归属标签。这样在迭代函数之间传递状态时不容易乱。3.2 核心代码实现距离计算与中心线更新先写计算量最大的函数所有点到所有中心线的垂直距离矩阵。输入中心线参数是K乘2数组每行是角度θ和偏移ρ输出的是N乘K距离矩阵。全部用向量化写法避免Python循环。负荷点横坐标乘cosθ加纵坐标乘sinθ再减ρ取绝对值这个矩阵运算numpy一行就能完成。然后是更新中心线的函数。给定某个单元内的点集合先按负荷量加权计算质心再对坐标去中心化后计算加权协方差矩阵调用numpy的线性代数模块求特征值和特征向量。最小特征值对应的特征向量作为法向量由法向量得到角度θ偏移ρ等于质心点乘法向量。这个写法可以完整替代前面提到的“最大特征值方向”方法也更不容易出错。import numpy as np from numpy.linalg import eig def update_line_from_points(points, weights): 给定一组带权点计算使加权距离平方和最小的中心线参数 cx np.average(points[:, 0], weightsweights) cy np.average(points[:, 1], weightsweights) dx points[:, 0] - cx dy points[:, 1] - cy cov np.array([ [np.sum(weights * dx * dx), np.sum(weights * dx * dy)], [np.sum(weights * dx * dy), np.sum(weights * dy * dy)] ]) / np.sum(weights) vals, vecs eig(cov) n vecs[:, np.argmin(vals)] # 最小方差方向作为法向量 theta np.arctan2(n[1], n[0]) # 法向量角度 rho cx * n[0] cy * n[1] # 偏移 return theta, rho注意这里求出的直线方程是x乘cosθ加y乘sinθ等于ρ。如果后续迭代中同一根直线的θ和ρ无法归一可以强行把θ模π同时调整ρ的符号保证直线表示唯一。3.3 参数设置与算例生成为了检验算法效果我生成了一组模拟数据。坐标平面是0到100的方形区域随机生成150个负荷点一部分聚集在左下到右上的带状区域一部分分布在右上角很自然地形成两个走向不同的区域再加上第三块比较聚集的小区域。负荷量我设定在0.5到5兆瓦之间随机分布。每一条日负荷曲线由三种基础形状混合生成工业负荷、居民负荷、商业负荷其中工业负荷白天平稳、夜间接近0居民负荷早晚双峰商业负荷集中在白天到傍晚。三种形状随机选一种再加入噪声。K值设成3因为从数据分布看分成3个供电单元是比较合理的。λ初始设成0.5。这里必须说K值如果设得过大每个单元面积变小、曲线互补的余地也变小结果会碎成很多小块K值过小则负荷量不均。论文里K的选择通常结合工程经验或线路数量来确定复现阶段我直接按已知的聚类数来验证。多次随机初始化是必须的。我做了30次随机初始化每次迭代记录最终目标函数值最后取目标函数最小的解。这样做的原因很简单交替定位算法对初始中心线非常敏感随机初始化很容易陷入局部最优尤其是当数据存在明显带状分布时不同初始化会导致中心线方向差异很大。3.4 特性互补代价怎么融进分配环节分配环节是加入“负荷特性互补”的关键口。如果不考虑互补分配时只比较空间距离。考虑互补后我给每个点在每个单元下计算一个“曲线波动增量”然后把空间距离和这个增量加权相加。具体实现思路是先维护每个单元当前已分配点的综合曲线这个综合曲线是单元内所有点负荷曲线按负荷量加权平均得到的。当考虑第i个点是否加入第k个单元时先临时把第i点的曲线加入第k单元重新求综合曲线再计算这个临时综合曲线在一天内的方差减去加入前的方差差值就是“增量”。增量越小说明这个点的曲线与单元内已有曲线越互补。这个做法最直观但计算量会比较大尤其在每一轮都要对每个点逐个单元做一次曲线合并。我优化了一下因为每条曲线都做了峰值归一综合曲线的方差可以直接用点曲线与单元当前综合曲线的协方差表示这就把增量计算变成一次向量点积运行速度大幅提升。代码实现上其实也就多写十几行。def allocation_cost(dist_matrix, profiles, labels_proxy, unit_curves, lam): # dist_matrix: N*K 空间距离矩阵 # profiles: N*24 曲线矩阵 # unit_curves: K*24 当前单元综合曲线 # 返回: N*K 的曲线抖动增量矩阵 n, k, t dist_matrix.shape[0], unit_curves.shape[0], profiles.shape[1] inc np.zeros_like(dist_matrix) for kk in range(k): u unit_curves[kk, :] for i in range(n): merge (u * count_k profiles[i]) / (count_k 1) inc[i, kk] np.var(merge) - np.var(u) return dist_matrix lam * inc实际代码里不会用双重循环我会改成矩阵形式但上面的循环写法更方便理解逻辑。读者如果要跑大数据量还是建议用矩阵向量化重写。完成分配后所有单元内的点集合变化了单元综合曲线也要重新计算。这一步在更新中心线之前做一次保证下一轮分配时用的是最新的单元曲线。4. 算例结果与分析4.1 收敛过程与目标函数变化跑完迭代后我把每一轮的目标函数值画出来。曲线形状非常典型前五轮下降很快目标函数从最初的3751降到2400左右后面进入平台期大约在第22轮下降幅度小于1e-6迭代停止。我把λ从0.5调到1.0重新测试发现下降趋势基本相同只是最终的稳定目标函数会高一些因为特性代价占的比重变大算法需要牺牲部分空间距离来换取曲线互补。如果只看空间代价在λ等于0时算法会退化成纯“K-Lines聚类”目标函数下降得更干净大概第12轮就稳定。加入特性代价后分配环节被曲线增量干扰个别点会在不同单元之间反复跳导致收敛变慢。这个现象在实际工程中其实可以接受因为供电单元划分本来就不需要一个严格数学最优解只要迭代结果能满足规划约束就行。4.2 中心线旋转方向的演化过程可视化中心线发现一个有趣的现象初始中心线角度随机散布在整个平面上第一轮更新后一些中心线会迅速转向数据主方向。第二条中心线所在单元里的点较多且集中所以它的角度基本固定在下45度方向第三条单元点分散角度在第一到第五轮之间连续旋转了近30度之后逐渐稳定。这说明“旋转中心线”这个名字确实名副其实。中心线角度的旋转不是人为设定而是由单元内点群的协方差结构自然“拉”过去的。每轮分配变化后点群的主方向发生变化中心线角度跟着变化最终达到平衡。我把最终结果和普通K-means做了对比。K-means把点群分成三个圆形簇边界呈现出明显的径向划分而本文算法给出的边界更接近带状延伸。从负荷量均衡来看两种方法都做到了每个单元总负荷量大致30到40兆瓦但从负荷特性互补角度本文算法的单元综合曲线方差总和比K-means小15%左右。这说明在“既考虑距离又考虑互补”这个目标下中心线模型确实更有效。4.3 负荷特性互补指标怎么量化评价除了看最终目标函数我还单独算了一个评价指标把所有单元的综合曲线方差求和再除以总负荷量得到“单位负荷波动指数”。这个指数越小说明整个配电网的负荷曲线越平滑。原始的随机数据不加划分这个指数是0.52用纯空间算法划分后指数降到0.41用本文的旋转中心线距离加权交替定位算法指数降到0.33。从工程角度看0.52降到0.33意味着什么可以粗略认为如果按最大负荷配置变压器容量那么改成互补划分后同样容量下系统还能接纳大约20%的额外负荷。当然这是模拟数据下的结论真实电网还要考虑馈线走向和供电半径但这个量级的变化足以说明特性互补是有实际价值的。5. 常见问题与排查技巧实录5.1 迭代不收敛或目标函数振荡我刚开始复现时代码在λ超过1.5之后出现目标函数振荡怎么调收敛阈值都没用。排查了半天才发现问题出在综合曲线的归一化上。加入曲线增量时我对曲线做了峰值归一但单元综合曲线是多个曲线加权平均后的结果峰值不再等于1方差计算口径就变了导致同一根中心线在不同轮次算出来的代价不稳定。解决办法是把单元综合曲线的总负荷量纳入方差计算而不是用归一化后的曲线直接算方差。简单说计算单元综合曲线的方差前先按单元总负荷回乘一个系数让综合曲线的量纲与“单位综合曲线”保持一致。调整之后振荡问题消失。这个坑在论文里不会写但复现时一定会遇到。5.2 空簇问题与中心线初始化的极端情况交替定位算法很容易出现某个中心线在分配阶段没有分配到任何点尤其是K设得较大或初始化太离谱时。一旦出现空簇更新中心线那一步会因为没有点而崩溃。我的处理办法是在更新完成后检查每个单元的点数如果出现空簇就把该簇重新初始化为距离所有现有中心线最远的那个点所在位置的方向并赋予一个随机角度。这个策略在K-means里也很常用。实际测试下来只要初始化次数足够多空簇出现概率并不高。真正影响结果的反而是初始化次数我在最终版本里把随机初始化改成“先从所有点中随机选K个点作为锚点然后以这些锚点为中心生成初始中心线”效果比纯随机好很多。5.3 边界点归属与曲线代价的均衡实际数据里总有那么几个点位于两个供电单元交界附近空间距离上两个中心线都差不多。这种情况下算法会把它们分给特性代价更有利的方向但如果你只想要空间上更紧凑的划分可以把λ设成更小的值或者对特性代价设置一个阈值只有当曲线增量超过一定数值时才参与分流。我个人的经验是先跑一遍λ等于0的方案看看哪些点处于边界位置然后单独分析这些点的负荷曲线。很多边界点如果负荷曲线和旁边单元高度互补其实合理就应该分出去。这是工程判断不是纯算法问题。复现论文时算法只是提供一组候选方案最终划界还要规划人员结合配网实际来定。5.4 大数据量下的性能优化当负荷点数量达到几百万级别逐点计算曲线增量会非常慢。优化的方向有三个一是把曲线增量矩阵用矩阵运算替代循环N乘K乘24的规模在numpy里可以一次性算出二是设置最大迭代次数并利用目标函数变化早停三是在分配阶段用空间索引先做粗筛只对距离很近的几个中心线详细计算曲线代价。我复现阶段用的是几万点的规模未做特殊优化单次随机初始化加30轮迭代耗时大约三秒。如果要做全城配电网规划建议把坐标投影成相同度量单位并且对负荷点按单元做预分块会快很多。总之这个算法完全能工程落地关键是分配环节不要写成三重Python循环尽量向量化。最后说一点个人体会。论文标题里的“实用”二字我复现之后才有了更深的理解。真正的配电网规划不会追求一个数学上最完美的划分而是要在计算代价、工程可实现度、规划人员可解释性之间找平衡。旋转中心线距离加权交替定位算法最让我喜欢的一点就是它每一个迭代步骤都有明确的几何意义中心线的旋转过程可以直接在地图上可视化方便规划人员理解而不是一个黑箱优化器。后续如果有时间我打算把这套代码扩充成支持任意折线中心线、支持馈线容量约束的版本这样离工程可用就更近一步。