☰
旋转中心线距离加权交替定位算法复现与供电单元划分实战
2026/10/2 18:35:29 网站建设 项目流程

最近折腾完一篇中压配电网供电单元划分论文的复现工作,算法名挺长,叫“旋转中心线距离加权交替定位算法”,放在配电网里解决的是负荷特性互补和供电单元划分问题。这类工作在实际规划中非常实用,但论文里的公式和流程往往描述得比较含蓄,复现起来会踩到不少坑。这篇文章我把整个复现过程、算法拆解和实操经验一次性讲清楚,想动手跑代码的可以直接照着走。

先说这个算法到底在干嘛。配电网规划中常要把一片区域内的负荷点划分成若干供电单元,每个单元未来对应一个电源点或变电站供电范围。传统做法就是聚类,比如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 交替迭代:划分与中心线更新

交替迭代的完整流程如下:

  1. 初始化每个单元的中心线角度θ_k,可以等间隔布置,也可以按负荷点分布密度设定,推荐后者。
  2. 固定所有θ_k,逐个负荷点计算隶属度,重分配归属。
  3. 固定分配结果,更新每个单元的边界角度。更新规则要让中心线向单元内所有负荷点的“角向重心”靠拢。具体来说,以O为极坐标原点,把所有属于单元k的负荷点转换到极角φ_i,用这些极角的圆统计均值(circular mean)作为新的边界朝向参考,再结合两侧相邻单元的中心线角度做约束平滑。
  4. 检查归属是否变化,如果变化量小于阈值或达到最大迭代次数,终止;否则回到第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 核心数据结构与初始化

我定义一个类来管理整个算法,数据容器用dataclass:

from 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], axis=0))[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[k+1] - angles[k] < 0.05: move = (0.05 - (angles[k+1] - angles[k])) / 2 angles[k] -= move angles[k+1] += move for k, unit in enumerate(units): unit.angle = angles[k]

这个更新方式并不是论文逐字对应的原始公式,因为不同论文对“中心线”的定义有差异,我这边是延续“过旋转中心旋转”的几何含义做的合理实现。你要是看的原论文用了不同的参数化方式,核心迭代骨架是一样的:先E步计算隶属度,再M步更新角度。

3.4 收敛判断与结果可视化

收敛判断我用两个指标:一是所有负荷点归属的变化数量,二是单元中心线角度的变化量。前者更直接,因为最终输出的是划分结果。我设置为:一轮迭代后,归属变化的点数少于总点数的0.5%,即认为收敛。

def has_converged(old_assign, new_assign, threshold=0.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-Means(k=4)和本文算法划分,结果如下表:

指标K-Means结果旋转中心线算法结果
最大峰谷差率0.720.58
单元同时率0.860.74
平均供电半径(km)1.942.11
单元负荷均衡系数1.321.11

K-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给不了你这种直观的边界线。

另外,如果你准备用这个算法写论文或者做工程方案,建议不要照搬我上面的简化实现。原论文里大概率有严谨的目标函数和约束表达式,你需要把那部分看透彻,再结合我的代码框架去理解每个变量到底在原式里对应什么位置。我上面给出的是一种工程可用的近似版本,真实场景下还需要根据你的负荷点数量、曲线采样点数、是否存在多电源候选点做适配。

最后再分享一个小技巧:跑任何迭代类算法之前,先用一个只包含三到五个负荷点的小样例把每一轮迭代的手算结果和代码结果对照一遍,确认你的距离加权方向和角度更新逻辑和论文一致。我当时就是因为角度正负号理解反了,导致前三天跑出来的结果都是镜像分布,还以为是算法缺陷,浪费了不少时间。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询