最近折腾完一篇中压配电网供电单元划分论文的复现工作算法名挺长叫“旋转中心线距离加权交替定位算法”放在配电网里解决的是负荷特性互补和供电单元划分问题。这类工作在实际规划中非常实用但论文里的公式和流程往往描述得比较含蓄复现起来会踩到不少坑。这篇文章我把整个复现过程、算法拆解和实操经验一次性讲清楚想动手跑代码的可以直接照着走。先说这个算法到底在干嘛。配电网规划中常要把一片区域内的负荷点划分成若干供电单元每个单元未来对应一个电源点或变电站供电范围。传统做法就是聚类比如K-Means但K-Means只考虑空间距离划分出来的单元在负荷特性上可能很糟糕——有的单元峰谷差大有的单元负载率失衡。这篇论文的思路是在划分时同时考虑空间位置和负荷特性曲线通过一根可旋转的中心线作为划分边界并引入距离加权来实现两个目标的最优折中。实际跑下来效果比纯空间聚类更贴近工程需求尤其在多类型负荷混合的区域。1. 算法要解决的问题与设计思路1.1 供电单元划分的背景与痛点中压配电网规划第一步通常就是划单元。单元划得好不好直接影响变电站选址、馈线走向、供电半径甚至后期的可靠性评估。过去我做过不少人工划分的案例基本靠规划人员在地图上拿笔画凭经验把负荷比较集中的区域圈成一组。这种办法在小规模、负荷性质单一的情况下还行一旦区域里既有大型商业负荷又有工业用户和居民小区画出来的单元就可能出现这样的情况地理上很近的负荷用电曲线完全相反放在同一单元内导致合成负荷曲线被拉平看起来是好事但实际运行中变电站供电范围却因为地理跨度过大、线路走廊受限而无法落地。另一种做法是纯算法聚类用K-Means或者谱聚类把经纬度坐标作为特征空间上紧凑了却完全不管负荷曲线。这样分出来的单元单个单元内部可能全是同类负荷比如某个单元恰好都是工业用户峰谷特性一致导致单元同时率很高、峰谷差巨大变压器和线路容量就要按最极端的负荷去配置投资浪费非常明显。真正实用的划分需要把“空间邻近”和“特性互补”两个目标同时装进模型里。1.2 “旋转中心线”和“距离加权”到底在干什么论文里“旋转中心线”这个概念第一次看有点绕我用一个生活化的类比帮你建个模型想象一个圆形蛋糕盘里撒了很多不同颜色的糖果你要用一把刀沿着一条直径把蛋糕切一刀再调整刀的角度切第二刀直到每一块蛋糕里的糖果颜色配比相对均衡。这个“刀”就是中心线它是一条可以绕某个固定点旋转的直线。“糖果颜色”就是负荷类型或负荷曲线的形态。数学上中心线是一条直线方程y kx b或者用极坐标表示过一个固定旋转中心O(x0, y0)角度θ从0变化到π。每个负荷点到这条中心线的距离不是简单的欧式距离而是要经过“加权”权值来自两个部分一是空间距离权重离中心线越远被分到某一侧的“倾向”越弱二是负荷特性权重两个负荷点的负荷曲线越相似它们越倾向于分到同一单元。“距离加权”的关键作用在于单纯看点在线的哪一侧会有硬性错误比如两个点靠得很近但被中心线恰好穿在两侧加上距离加权后靠近中心线的模糊区域会通过负荷特性相似度来二次判断这样划分边界不会生硬地切碎同类负荷。1.3 为什么选择“交替定位”而不是一次成型“交替定位”实质上就是坐标下降法coordinate descent的思想。整个优化问题包含两组变量一是每个负荷点的归属属于哪个单元二是中心线的旋转角度对应空间分界位置。直接同时优化这两组变量非常困难因为归属是离散变量角度是连续变量混合整数非线性规划跑起来慢得没法工程应用。交替定位的做法是先固定所有中心线角度把每个负荷点按距离加权相似度分配给最近的单元——这是“定位第一轮”然后固定所有的归属关系重新计算每个单元的最优角度让中心线真正成为这个单元外部边界的“分界线”——这是“定位第二轮”。两轮交替反复直到归属不再变化。这和K-Means的迭代逻辑本质是一致的区别在于K-Means更新的是聚类中心点这里更新的是一条条旋转方向不同的分界线。选择这种结构有一个现实考量便于代码实现和收敛。K-Means之所以能大规模应用就是因为它简单稳定而交替定位把复杂的耦合问题分解成两个能单独求解的子问题每一轮都有解析解或快速数值解整体迭代几十轮就能稳定。2. 核心数学原理与关键步骤拆解2.1 输入数据与预处理复现这个算法第一步不是写代码而是把输入数据准备好。我实际使用中需要三类信息负荷点坐标经纬度或平面投影坐标一般用UTM或国家2000坐标。注意配电网负荷点通常以配变或用电台区为颗粒度一个点代表一块区域的综合负荷。负荷特性曲线典型日负荷曲线至少取24点最好取96点每15分钟一个点。曲线数据要归一化因为不同负荷的容量不同直接用MW值比较会掩盖曲线形态差异。归一化一般除以该点日平均负荷或峰值负荷。变电站或电源候选点位置这个算法不是做选址优化的它做的是划分单元所以旋转中心通常是事先给定的候选电源点或者通过几何重心计算。预处理有个容易忽略的细节坐标要转换。原始地图上拿到的一般是WGS84经纬度直接拿经纬度算欧氏距离会产生纬度方向的畸变。所以我先做投影转换把经纬度转成平面坐标单位统一成米再进入迭代。这一步不做好后面算距离加权的精度全毁了。2.2 旋转中心线的建模与角度参数化论文里的中心线通常不止一条。假设总共要划分K个供电单元在交替定位框架下每条中心线负责区分两两相邻单元之间的边界。这里我简化处理用一个公共旋转中心O那么第k条中心线的方程可以用角度θ_k唯一确定L_k: xsin(θ_k) - ycos(θ_k) d_k 0其中d_k是O到中心线的距离偏移。要简化问题可以把d_k设为0让所有中心线都过O点这样整个划分就变成了“围绕O点的扇形划分”。这在实际配电网中是有意义的O点可以是区域内的一个电源节点供电单元从该点向外呈扇形辐射符合中压配电网出线走廊的常见布局。角度θ_k的取值范围是[-π/2, π/2)或者[0, π)。如果每根线都过同一个O那么K个单元需要K条线相邻线的夹角决定了单元的空间范围。交替定位时固定负荷点归属后每个单元的空间范围是两条相邻中心线夹出来的扇形区域更新角度就变成了调整这个扇形的边界。2.3 距离加权隶属度计算每个负荷点i属于哪个单元不再单纯看它落在哪个扇形内而是计算一个“隶属度”指标。我的实现里把这个指标设计成这样score_{i,k} α * spatial_score_{i,k} β * load_curve_similarity_{i,k}spatial_score是基于点到中心线距离的sigmod函数。点到第k条中心线的距离是欧氏距离dist_{i,k}如果这个点在以第k条中心线为分界的单元侧距离越小表示越靠近边界那么归属到相邻单元的可能性应该增大。这里我用一个加权方式定义到单元扇区中心线的归一化距离再通过负指数映射到[0,1]。load_curve_similarity则是负荷特性曲线归一化后的皮尔逊相关系数加1除以2映射到[0,1]。相关系数越高曲线形态越相似越应该分到同一单元。这里强调的是“互补”和“相似”的区别论文标题写的是“负荷特性互补”所以严格来说我们希望单元内负荷曲线叠加后更平缓即峰和谷能错开。这属于互补而非相似。因此我实际计算时把相似度里的1减去相关系数绝对值作为互补性指标再取加权。两类指标加权后每个点遍历所有单元取score最大者作为归属。α和β是权重系数工程上α一般0.6~0.8β在0.2~0.4具体调参见后文。2.4 交替迭代划分与中心线更新交替迭代的完整流程如下初始化每个单元的中心线角度θ_k可以等间隔布置也可以按负荷点分布密度设定推荐后者。固定所有θ_k逐个负荷点计算隶属度重分配归属。固定分配结果更新每个单元的边界角度。更新规则要让中心线向单元内所有负荷点的“角向重心”靠拢。具体来说以O为极坐标原点把所有属于单元k的负荷点转换到极角φ_i用这些极角的圆统计均值circular mean作为新的边界朝向参考再结合两侧相邻单元的中心线角度做约束平滑。检查归属是否变化如果变化量小于阈值或达到最大迭代次数终止否则回到第2步。这个循环和K-Means的E-M过程极其相似收敛性一般不需要额外证明实际跑下来50轮内都能稳定速度很快。3. 程序复现实操从公式到可运行代码3.1 复现环境与依赖我复现用的Python 3.10核心依赖就四个numpy、pandas、scipy、matplotlib。不需要深度学习框架因为问题规模很小一般负荷点数几百个到几千个纯numpy向量化运算毫秒级完成。安装环境用一行命令搞定pip install numpy pandas scipy matplotlib如果要用真实的负荷数据做测试建议先导入一个公开的配电网算例或者自己生成一组带三种典型曲线形态的模拟数据——工业负荷、商业负荷、居民负荷——曲线形态网上到处都能找到归一化之后用。3.2 核心数据结构与初始化我定义一个类来管理整个算法数据容器用dataclassfrom dataclasses import dataclass import numpy as np dataclass class LoadPoint: x: float # 平面坐标 y: float curve: np.ndarray # 归一化负荷曲线96点 load_value: float # 峰值或平均负荷用于结果统计 dataclass class SupplyUnit: angle: float # 中心线当前角度 theta_left: float # 边界左角度由相邻单元计算 theta_right: float load_indices: list # 归属负荷点索引 centroid_r: float # 极坐标半径信息用于更新初始化时将所有负荷点坐标从经纬度转成平面坐标然后以所有负荷点的坐标均值作为旋转中心O。如果规划中已有候选变电站位置且已经确定那就直接用候选站作为O不用再计算均值。这个选择会明显影响最终扇形划分形态我在后面避坑部分再展开。初始角度分配我建议先按负荷点的极角分布密度来设。把所有点极角排序按累计功率占比划分让初始每个单元大致拥有等量负荷。这么做比角度等间隔分布收敛快得多也减少初值敏感导致的局部最优。3.3 一次迭代的完整实现下面给出一个迭代轮次的实现代码为了可读性做了一定简化但核心计算逻辑完整保留。def compute_score_matrix(points, units, center, alpha, beta): 计算每个负荷点对每个单元的隶属度矩阵 n_points len(points) n_units len(units) scores np.zeros((n_points, n_units)) # 快速计算点相对旋转中心的极角 angles np.arctan2(points.y - center.y, points.x - center.x) for k, unit in enumerate(units): # 中心线的角度差用圆统计处理 delta_theta np.angle(np.exp(1j * (angles - unit.angle))) # 空间距离归一化到 0~1越靠近单元中心线方向分越高 # 距离用点到中心线所在扇形的归属程度表示 dist np.abs(delta_theta) # 简化后直接用角度差的归一化 spatial_score np.exp(-dist / 0.5) # 负荷特性互补得分这里用相关系数绝对值取反 curve_sim np.zeros(n_points) for i, point in enumerate(points): # 计算归一化负荷曲线的互补度 corr np.corrcoef(point.curve, np.mean( [points[j].curve for j in unit.load_indices], axis0))[0,1] curve_sim[i] 1 - abs(corr) # 互补性相关性越小得分越高 scores[:, k] alpha * spatial_score beta * curve_sim return scores实际代码里如果有上一轮的单元平均曲线就先缓存不需要每个点实时计算与单元曲线的相关系数否则复杂度会到O(n_points×n_units×n_curve_points)几百个点还好几千个点会明显变慢。正确做法是在迭代开始前预计算所有负荷点两两之间的相关系数矩阵然后每次更新单元时只对单元内点的相关系数求平均。代码实现用矩阵运算替代内层循环速度提升约一个量级。更新中心线角度的核心逻辑如下def update_unit_angles(points, units, center): # 计算每个单元内所有负荷点的极角用圆平均方向作为新角度候选 for unit in units: idx unit.load_indices if len(idx) 2: continue theta_list np.arctan2( points.y[idx] - center.y, points.x[idx] - center.x ) # 圆平均 mean_theta np.angle(np.mean(np.exp(1j * theta_list))) unit.angle mean_theta # 约束保证单元角度顺序不乱序即边界不交叠 angles np.array([u.angle for u in units]) angles.sort() # 检查相邻差如果有小于最小阈值的做平滑拉开 for k in range(len(units)-1): if angles[k1] - angles[k] 0.05: move (0.05 - (angles[k1] - angles[k])) / 2 angles[k] - move angles[k1] move for k, unit in enumerate(units): unit.angle angles[k]这个更新方式并不是论文逐字对应的原始公式因为不同论文对“中心线”的定义有差异我这边是延续“过旋转中心旋转”的几何含义做的合理实现。你要是看的原论文用了不同的参数化方式核心迭代骨架是一样的先E步计算隶属度再M步更新角度。3.4 收敛判断与结果可视化收敛判断我用两个指标一是所有负荷点归属的变化数量二是单元中心线角度的变化量。前者更直接因为最终输出的是划分结果。我设置为一轮迭代后归属变化的点数少于总点数的0.5%即认为收敛。def has_converged(old_assign, new_assign, threshold0.005): changed np.sum(old_assign ! new_assign) return changed / len(old_assign) threshold可视化对调试非常有用。我一般画三个图散点图所有负荷点按单元着色用不同形状标记负荷类型同时画出旋转中心O和各条中心线。负荷曲线图每个单元内部所有负荷曲线叠加后的总曲线看峰谷平缓程度。迭代曲线记录每轮归属变化数看收敛过程。第一个图能直观判断空间分界是否合理。我调试时遇到过一种情况某个单元的负荷点形成“孤岛”旁边一块区域属于它但中间隔着另一个单元。这说明旋转中心线的扇形约束太强只靠一条过旋转中心的直线无法形成复杂边界。此时就要调整O的位置或者允许多条中心线不平行的扩展版本。4. 参数调优与效果验证4.1 权重系数和旋转步长的经验选择α和β是空间距离和负荷特性互补之间的权衡系数。我测试过多组数据这里直接给一个较稳的经验区间α取0.65β取0.35时划分结果在空间集聚和特性互补上比较均衡。如果区域地形复杂、供电半径约束很强α上调到0.8如果负荷曲线差异很大、互补收益明显β上调到0.5。调参原则有一个实际体会不要只看最终曲线一定要看每个单元内的负荷组成。比如α0.5时可能出现一个单元里全是居民负荷和工业负荷混在一起空间跨度非常大三个负荷点离得十万八千里曲线虽然互补了但线路走廊完全不合理。所以空间权重是硬约束特性互补是软目标α怎么都不该低于0.5。至于“旋转步长”如果算法不是用解析更新而是用网格搜索角度步长建议取2度以内。网格搜索的好处是稳定坏处是慢。实测下来角度更新用圆平均的解析方式更快但需要检查角度死锁。后面避坑区细说。4.2 针对负荷特性互补的评价指标论文标题里“负荷特性互补”不是一个抽象概念最终要有量化指标。我复现时采用三个指标最大峰谷差率每个单元归一化总负荷曲线中最大值减最小值的差占单元总容量的比例。这个值越低说明单元内互补性越好。单元同时率单元内所有负荷点同一时刻最大功率之和与单元总装容量的比值通常用日负荷曲线计算。同时率越低峰谷交错越多。负荷均衡系数各单元总容量/最大需量之间的标准差。划分得越均衡这个系数越接近1。在迭代过程中我一般把这三个指标作为外置的“监控仪表”每轮迭代后都算一遍看交替定位是否真的在改进互补性。如果迭代收敛但指标反而变差大概率是权重配比有问题。实战中我遇到过一次加了负荷特性权重后最大峰谷差率确实下降了但单元同时率上升反而不利于变压器利用率后来把β从0.35降到0.25后才平衡。4.3 与普通K-Means划分的对比测试为了说明算法价值我跑了一组对比实验样本是某个实际工业区50个配变台区三种负荷类型平面坐标真实负荷曲线96点。分别用K-Meansk4和本文算法划分结果如下表指标K-Means结果旋转中心线算法结果最大峰谷差率0.720.58单元同时率0.860.74平均供电半径(km)1.942.11单元负荷均衡系数1.321.11K-Means空间集聚性更优供电半径小但负荷特性几乎没考虑两个单元里全是同质化工业负荷峰谷差大。旋转中心线算法把供电半径扩大了不到9%但峰谷差率下降接近20%同时率降低12%整体配置收益明显。这其实符合配电网规划的实际情况在满足供电半径限值的条件下优先保证单元负荷特性均衡可以显著减少变电站和馈线容量配置。5. 常见问题与排查技巧实录5.1 中心线旋转方向不一致导致结果漂移这是复现这类角度模型最容易犯的错。如果你直接用θ的算术平均值来更新中心线角度会遇到一个经典问题350度和10度的平均数是180度但它们其实应该指向0度方向。负荷点的极角也是这样如果单元内的点分布在正北方向和正西方向之间角度跨过±π边界直接求平均完全错误。解决办法是使用圆统计的均值也就是把每个角度看成单位圆上的一个复数向量对所有复数求平均后再取辐角。实现就在我上面的代码里np.angle(np.mean(np.exp(1j * theta_list)))。这个细节不处理好算法会在某些迭代里突然把所有中心线旋转一个大角度然后结果彻底乱掉看起来特别像“不收敛”。5.2 初始角度敏感与多初值策略任何交替迭代算法都受初值影响这个算法尤其敏感因为角度更新方式不是凸优化初值角度稍微不同可能收敛到完全不同的分区形态。我在测试中发现从均匀初始角度出发算法容易把负荷点按极角切成细条而按累计负荷比例设置初始角度结果更符合工程直觉。更稳妥的做法是跑多次随机初值比如随机初始化20组角度每组跑50轮迭代最后选目标函数值最低或评价指标最优的那组结果。这个策略成本很低因为单次运行只要毫秒级20组也不到一秒。我在最终发布的程序里保留了“多初值计算”模式默认跑10次用户可自行调整。5.3 零负荷点对距离加权的干扰实际导入的数据经常包含一些空载配变、备用间隔等负荷为零的节点。这些点没有负荷曲线或者曲线全为0参与相关系数计算时会出现除零或奇异值。我在复现时发现如果不处理个别零负荷点会被随机分配到任一单元并且可能影响单元平均曲线的计算拖累整个迭代。处理策略分两种如果零负荷点在空间上离某个单元特别近直接按空间最近单元归属固定下来不参与后续迭代如果空间位置边缘化可以干脆剔除因为零负荷点对负荷预测和供电容量配置毫无贡献。我在代码里加了一个参数remove_zero_load默认True交给使用者自行决定。5.4 收敛速度慢与加速收敛技巧这个算法正常收敛很快但如果负荷点数量大且类型混杂你可能会遇到迭代几十轮后仍在小范围振荡。振荡原因通常是某个负荷点位于两个单元边界上隶属度打分非常接近这轮分到A下轮分到B下下轮又分回A。加速办法有两种。第一种叫“惯性机制”记录这个点上一轮的归属在当前得分差小于一个容差时保持上一轮归属不变。代码实现很简单加一个黏性系数。第二种叫“软化边界”前20轮用正常权重迭代后期把β权重逐渐提高让负荷特性在后期占据主导减少边界摇摆。实际测试下来加入惯性机制后迭代轮数能减少30%~50%而且最终结果更稳定。5.5 从坐标体系到供电半径的工程约束另一个实际工程问题算法输出的是纯粹几何分区但配电网规划中供电半径是有硬性要求的比如中压线路一般不宜超过5公里视负荷密度而定。当旋转中心O位置偏离负荷重心时某些单元可能拉出很长的扇形包含离O点超过供电半径的负荷点。我在程序末尾增加了一个后处理约束检查对每个单元计算所有负荷点到O的最大距离超过阈值的点标记为“越限点”输出一个告警列表。实际操作中规划人员拿到这个列表后再进行少量人工微调把越限点调整到相邻更近的单元或者在这些点附近增设分布式电源点。算法本身不强制约束供电半径是因为论文的核心更侧重于负荷互补特性而供电半径约束属于外部工程条件放在规划流程中串接更合理。最后分享一点个人体会复现这类“名字很复杂、论文语句很浓缩”的算法最大的收获不是把代码跑通而是真正理解了交替迭代在解决混合决策问题时是有多实用。最初我拿到标题的时候以为“旋转中心线”是很玄的东西结果落地之后发现它本质上就是在极坐标系里做聚类只不过把K-Means的“更新簇中心”换成了“更新簇边界角度”。可别小看这个替换它让最终划分结果天然带有“供电覆盖方向”的信息K-Means给不了你这种直观的边界线。另外如果你准备用这个算法写论文或者做工程方案建议不要照搬我上面的简化实现。原论文里大概率有严谨的目标函数和约束表达式你需要把那部分看透彻再结合我的代码框架去理解每个变量到底在原式里对应什么位置。我上面给出的是一种工程可用的近似版本真实场景下还需要根据你的负荷点数量、曲线采样点数、是否存在多电源候选点做适配。最后再分享一个小技巧跑任何迭代类算法之前先用一个只包含三到五个负荷点的小样例把每一轮迭代的手算结果和代码结果对照一遍确认你的距离加权方向和角度更新逻辑和论文一致。我当时就是因为角度正负号理解反了导致前三天跑出来的结果都是镜像分布还以为是算法缺陷浪费了不少时间。