电动汽车、主动配电网与电力系统规划下的Voronoi图围捕算法Matlab复现
说实话,第一次看到这个题目的人大概率会被绕晕——电动汽车、主动配电网、电力系统规划、Voronoi图、围捕算法,五个词单拎出来都认识,拼在一起就不知道是个什么东西了。
简单说,这是我前段时间复现的一个偏研究性质的项目,核心思路是:在配电网规划问题里,用Voronoi图做空间分区,再用围捕算法做优化搜索,解决电动汽车充电负荷接入后,主动配电网的站址选址、供电分区和资源调度协同问题。听起来很高端,落地其实就是一个Matlab仿真工程,全套东西跑通之后,我最大的感受是:这个组合思路确实有意思,但网上资料零散,很多细节要靠自己踩坑试出来。
这篇文章我尽量把整套逻辑、数学建模、代码实现、参数设置和常见的坑都讲透,适合正在做配电网规划、电动汽车充电设施布局、或者对Voronoi这类空间分区算法在工程优化里怎么用感兴趣的朋友。哪怕你是刚接触这块的初学者,只要会用一点Matlab,照着做也能跑出自己的结果。
1. 项目概述与整体思路拆解
1.1 这个项目到底在解决什么问题
先说背景。传统配电网规划一般考虑的是负荷增长和网架扩展,核心是“在哪建变电站、怎么接线最省钱”。但电动汽车大规模接入之后,问题性质变了——充电负荷不再是简单的“预测数值”,它带有强烈的空间分散性和时间随机性。
比如某个区域白天写字楼附近充电需求大,晚上居民区附近需求大;同一个站址,上午可能排队、下午可能闲置。这种时空错配如果不在规划阶段处理,后面运行阶段怎么调度都费劲。
主动配电网的概念就是在这种背景下提出的,核心是让配电网从“被动接受负荷”变成“主动引导资源”。而规划层面就要回答:充电站建在哪、每个站的供电区域怎么划、分布式电源怎么配合、网架怎么改。这些问题本质上是一个空间优化问题,Voronoi图在这里的价值就是自动完成区域划分,把一个连续地理空间切割成若干个多边形,每个多边形对应一个设施的供电范围或服务范围。
1.2 为什么要把Voronoi图和围捕算法组合在一起
这是整个项目最容易被问的问题。我当时也疑惑:Voronoi图不是几何划分工具吗?围捕算法不是机器人路径规划里的东西吗?这两者怎么凑到一起?
答案在于优化问题的结构。配电网规划是个典型的混合整数非线性规划,变量包括:站址坐标(连续)、站容(离散)、馈线拓扑(0-1)、分布式电源位置(0-1),目标函数包含投资成本、网损、可靠性、碳排放等多个维度。直接求解基本不可能在可接受时间内完成。
而分层或者分解思路是常用的降维方式。第一层用遗传算法、粒子群这类元启发式搜索站址;第二层在给定站址下,用 Voronoi 划分服务区域、确定每个站服务的负荷;第三层做潮流计算评估方案优劣。这是比较经典的分层规划套路。
但单纯用遗传算法搜站址有个问题:收敛慢、容易陷入局部最优,尤其当区域里有多个充电需求热点时,站址很容易扎堆,Voronoi划分出来就会出现“一个站忙死、一个站闲死”的方案。
围捕算法这时候就派上用场了。它的思想来自多智能体协同围捕:多个捕食者从不同方向逼近猎物,彼此保持角度间隔,逐步收缩包围圈。如果把充电站当成捕食者、充电需求热点当成猎物,那“围捕”的过程天然保证了站址的空间分散性,同时动态调整站址朝向需求中心移动。这正好解决了我刚才说的扎堆问题。
所以这个组合不是生拼硬凑,而是用Voronoi的几何划分能力做空间解耦,用围捕策略的运动学特性做站址搜索,两者在数学上形成互补。
1.3 完整的技术链路预览
我在复现时把整个流程拆成了七个模块,这里先给一个全景图:
- 基础数据生成:区域道路/负荷坐标、电动汽车充电需求(蒙特卡洛抽样)、分布式光伏出力曲线。
- 场景建模:把研究区域离散成负荷点,每个负荷点带充电需求权重。
- Voronoi分区:给定一组站址,生成Voronoi图,把负荷点分配到对应的站点服务区。
- 围捕算法:把当前站址集合看作捕食者群体,以负荷重心为猎物,迭代更新站址位置。
- 目标函数评估:对每个分区方案计算年综合费用(投资+网损+运行+碳排放),调用潮流计算。
- 外层寻优:围捕迭代结束后,再叠加一次粒子群微调,避免收敛到围观意义上的局部最优。
- 可视化与结果分析:画出Voronoi图、围捕轨迹、收敛曲线、不同方案对比图。
整个过程用Matlab R2021b以上版本可以完整实现,不需要额外工具箱,核心几何计算用polyshape和自写的Voronoi裁剪函数,优化求解只用基础的矩阵运算。
2. 算法原理与核心细节解析
2.1 Voronoi图原理:从“势力范围”到供电分区
Voronoi图说起来很简单:平面上给定N个种子点,把整个平面划分成N个区域,每个区域包含距离该种子点最近的所有点。用生活化的语言讲,就是每个种子点画一个“势力范围”,范围内的任何位置,到这个种子点的距离都比到其他种子点近。
数学表达式长这样:
区域 V_i = { x ∈ R^2 | dist(x, p_i) ≤ dist(x, p_j), ∀ j ≠ i }
这里的 dist 可以是欧氏距离,也可以根据实际场景换成路网距离、加权距离。
在配电网规划里用Voronoi图,我的理解是它天然做了三件事:
- 空间解耦:把一个大区域拆成若干子区,每个子区只跟一个供电设施挂钩,全局优化问题就变成了多个独立子问题,计算复杂度显著下降。
- 边界确定:相邻站的供电边界自动确定,不需要手动划分,且边界是直线段。这比人工画边界要客观得多。
- 负荷就近分配:如果负荷到站点的距离就是电气距离的近似,那Voronoi分区本质上是在最小化“负荷到服务设施的加权距离总和”,对应到规划里就是降低低压侧线损。
不过直接调Matlab自带的voronoi()函数有个大坑:它会生成无界区域。四周边上的种子点对应的Voronoi多边形边界跑到无穷远处去了。这在数学上没问题,在工程里就不行——研究区域是有边界的。所以必须做有界Voronoi裁剪,把无界多边形限制在研究矩形范围内。后面我会详细说这个怎么处理。
2.2 围捕算法的核心逻辑:需求热点牵引下的站址更新
围捕算法的原始场景是我方多个智能体捕捉一个移动目标,每个捕食者运动时考虑两个因素:朝向目标的引力和与同伴之间的斥力。引力让捕食者靠近目标,斥力让捕食者保持相同的角度间隔,避免扎堆。
在配电网规划问题里,我把这个机制做了一个映射:
- 猎物 = 分区内负荷重心(或充电需求热点位置)
- 捕食者 = 充电站/变电站候选站址
- 引力项 = 站址向本区负荷重心移动,保证服务效率
- 斥力项 = 站址之间保持距离,避免供电范围叠加重叠
- 包围收缩 = 每次迭代后需求热点变化,各站重新计算重心、重新逼近
具体到每个站的更新公式,我用的是改进的人工势场形式:
新位置 = 当前位置 + α × (负荷重心 - 当前位置) - β × Σ(与相邻站的距离小于安全阈值时的排斥力方向)
这里 α 是引力步长系数,β 是斥力步长系数。α 太大站址会震荡,太小收敛慢;β 太大站址会过度分散,导致部分负荷距离过远。
一个比较关键的细节是:围捕迭代不是一步到位的,而是分区-计算重心-移动-重新分区循环。因为站址动了,Voronoi区域就变了,区域变了重心也跟着变。这个循环实际上是坐标下降法的几何版本,收敛性和稳定性在多次实验里表现不错。
2.3 为什么选Matlab而不是Python
用Matlab做这个项目,最重要的原因是配电网潮流计算生态比较成熟。我自己比较熟悉的是用稀疏矩阵手写牛顿-拉夫逊法潮流,Matlab的矩阵运算效率天然占优。另外polyshape对象做多边形布尔运算非常方便,裁剪Voronoi图、计算区域面积、判断负荷点是否落在某个多边形内都是几行代码的事。
Python也有scipy.spatial.Voronoi和shapely可以做类似的事,但我个人觉得Matlab在做这个场景时调试体验更好:变量窗口直接看数组结构,画图交互方便,跑优化循环的时候不用频繁print进度。
如果你最终要用Python部署,我的建议也是先用Matlab把算法逻辑验证清楚,再翻译过去,不要一上来就在Python里纠结几何库的边界行为。
3. 系统建模与数学描述
3.1 研究区域离散化与负荷模型
我先说数据层怎么建模。为了保证仿真可复现,研究区域我用了某区域简化配电网场景,一个10×10km的矩形区域,内部随机生成300个负荷节点,每个节点带坐标 (x, y) 和负荷类型标签(居民/商业/工业)。
电动汽车充电负荷不能简单当成固定值,我用蒙特卡洛抽样生成:
- 每辆电动汽车的起始充电时间服从早晚双峰分布(早上9点一个峰、晚上19点一个峰)
- 充电功率取7kW(慢充)和60kW(快充)两种
- 每天的充电需求量按行驶里程抽样,再折算成充电时长
叠加上这个充电负荷后,每个节点的综合负荷 = 基础负荷 + 电动汽车充电负荷期望值。这一步是后续Voronoi分区的权重输入。
这里有一个容易犯的错:直接用平均负荷做规划会掩盖峰值问题。我建议做三个场景:平日、节假日、极端高温日,把三个场景下的最大负荷作为规划负荷。这样规划出来的站容有裕度,不会出现实际运行中过载。
分布式光伏出力我也做了简化建模,采用Beta分布描述辐照度的随机性,再折算成功率输出。光伏接入的位置主要选在负荷密度较高的Voronoi区域内,这样光伏就地消纳比例高、反送功率小。
3.2 目标函数与约束条件
目标函数是年综合费用最小,包含四个部分:
年综合费用 = 投资等年值 + 运行维护费用 + 网损费用 + 碳排放成本
投资等年值是把充电站建设投资和配电线路建设投资按折现率换算成等年值。这里有个简单公式:
等年值 = 总投资 × r(1+r)^n / ((1+r)^n - 1)
其中 r 是折现率,n 是设备寿命年限。我用的折现率是8%,充电站寿命按15年算,线路按20年算。
网损费用是每次潮流计算之后统计系统总网损,再乘以电价折算成年费用。碳排放成本根据网损对应的等值碳排放量和碳价计算。
约束条件包括:
- 节点电压约束:0.95 ≤ U_i ≤ 1.05 pu
- 支路容量约束:S_ij ≤ S_ij^max
- 充电站服务能力约束:站内充电桩数量 ≥ 本Voronoi区内充电需求的峰值排队量
- 供电半径约束:任意负荷点到其所属站点的距离 ≤ 允许最大供电半径(本项目取3km)
- 光伏接入容量约束:不超过该节点变压器容量的30%
这些约束在Matlab里通过潮流计算返回结果、再附加惩罚函数的方式处理。违反约束的方案直接给一个很大的惩罚费用,让寻优算法自动避开。
3.3 双层优化模型的结构
整个优化整体是双层的:
外层是围捕算法+粒子群微调的混合寻优,决策变量是站址坐标集合 {(x_1, y_1), (x_2, y_2), ..., (x_K, y_K)}。
内层是Voronoi分区+潮流评估,给定站址后,把所有负荷点分配到最近的站点,在每个分区的边界上设定虚拟联络开关,按照辐射状网络跑潮流,返回目标函数值和约束违反量。
之所以做双层而不是单层一步到位,是因为原始问题里站址和供电范围是耦合的——站址变了,最优供电范围就变了;供电范围变了,潮流结果也变了。如果放在一起求解,变量维度太高,非线性太强,基本算不动。双层结构的好处是每次内层计算都能得到当前站址下的“最优分区”,外层只负责探索不同的站址组合,这个思路在工程优化里非常常见。
4. Matlab复现完整流程
4.1 环境准备与基础数据生成
我用的版本是Matlab R2021b,Windows 11环境运行,16GB内存跑300个节点、双层优化100轮、每轮10次潮流计算,耗时大约15分钟。如果你机器性能一般,建议把节点数降到150个,逻辑不变,跑起来更快。
第一步生成研究区域和负荷点:
% 研究区域范围 region_x = [0, 10]; % km region_y = [0, 10]; % km % 随机生成负荷节点 rng(42); % 固定随机种子,保证可复现 n_nodes = 300; load_nodes = rand(n_nodes, 2); load_nodes(:, 1) = load_nodes(:, 1) * (region_x(2) - region_x(1)) + region_x(1); load_nodes(:, 2) = load_nodes(:, 2) * (region_y(2) - region_y(1)) + region_y(1);注意这里固定了随机种子,否则每次跑随机数据不一样,对比方案没有意义。如果你做敏感性分析,可以改成循环跑多个随机种子,统计平均值。
第二步生成电动汽车充电需求,我用的是双峰正态分布抽样:
% 充电起始时间双峰分布 num_ev = 5000; peak_morning = 9; % 上午9点 peak_evening = 19; % 晚上19点 sigma = 1.5; t_start = zeros(num_ev, 1); for i = 1:num_ev if rand < 0.45 t_start(i) = peak_morning + sigma * randn; else t_start(i) = peak_evening + sigma * randn; end end % 剔除24点以后的数据,超出部分往前翻 t_start(t_start > 24) = t_start(t_start > 24) - 24;然后算每个负荷节点的充电负荷期望值,需要统计每个节点附近的充电需求。我用的方法是把每辆EV按坐标分配到最近的负荷节点(欧氏距离),再累加充电功率。
这个思路比较简单,但实测下来够用。如果你要严谨,应该用路网距离或者考虑充电站容量对用户选择的影响,不过那属于另一个层面的研究课题了。
4.2 Voronoi有界裁剪:绕开Matlab自带函数的坑
Matlab自带voronoi(x, y)只能画图,要拿到每个区域的顶点坐标,得用voronoin。但前面说了,直接生成的无界区域没法直接用。
我推荐的处理方法是:先用标准Voronoi拿连接关系,再用矩形边界裁剪每个多边形。核心代码可以这样写:
function polys = bounded_voronoi(sites, bounds) % sites: Kx2 站址坐标 % bounds: [xmin xmax ymin ymax] % 返回每个区域对应的polyshape对象 K = size(sites, 1); [V, C] = voronoin(sites); xmin = bounds(1); xmax = bounds(2); ymin = bounds(3); ymax = bounds(4); polys = polyshape.empty(K, 0); for k = 1:K % 取出第k个区域的顶点 region_verts = V(C{k}, :); % 去掉Inf顶点,这些顶点对应无界区域 infinite_idx = any(isinf(region_verts), 2); finite_verts = region_verts(~infinite_idx, :); % 如果有Inf顶点,需要构造一个足够大的裁剪框 if any(infinite_idx) % 扩展研究区域,确保裁剪后覆盖完整 pad = 10; ext_bounds = [xmin-pad xmax+pad ymin-pad ymax+pad]; ext_box = polyshape([ext_bounds(1) ext_bounds(1) ext_bounds(2) ext_bounds(2)], ... [ext_bounds(3) ext_bounds(4) ext_bounds(4) ext_bounds(3)]); % 对ext_box再进行一次Voronoi裁剪 % 简便做法:直接把有限顶点加上扩展框的四个角 ext_corners = [xmin-pad, ymin-pad; xmin-pad, ymax+pad; ... xmax+pad, ymax+pad; xmax+pad, ymin-pad]; all_verts = [finite_verts; ext_corners]; else all_verts = finite_verts; end % 创建多边形并裁剪到研究区域 shp = polyshape(all_verts(:, 1), all_verts(:, 2)); region_box = polyshape([xmin xmin xmax xmax], [ymin ymax ymax ymin]); polys(k) = intersect(shp, region_box); end end这个函数是整套代码的地基,也是坑最多的部分。我踩过的问题包括:polyshape在顶点共线时报警、两个多边形相交精度异常、边界上的站址产生退化多边形。后来我的处理方式是:
- 站址生成时加一个最小间距约束,保证任意两个站址距离不小于0.3km,避免退化。
polyshape操作后调用polys(k) = simplify(polys(k)),把退化部分清理掉。- 如果某个区域面积为0或顶点数少于3,直接跳过或回退到上一次迭代的站址。
4.3 围捕策略的核心实现
围捕迭代是算法的心脏。我实现的是一个简化版的多捕食者协同,每一步做三件事:计算各分区负荷重心 → 更新站址(引力+斥力) → 重新分区。
计算负荷重心时,把节点负荷作为权重:
function centroid = load_centroid(nodes_in_region, loads_in_region) total_load = sum(loads_in_region); if total_load <= 0 centroid = mean(nodes_in_region, 1); else centroid = sum(nodes_in_region .* loads_in_region, 1) / total_load; end end更新站址的核心代码:
% 参数设置 alpha = 0.35; % 引力系数 beta = 0.15; % 斥力系数 safe_dist = 2.5; % 站间安全距离(km) max_iter = 30; % 初始化站址(用K-means跑一次得到初始点,避免初始站址扎堆) [~, init_idx] = maxk(rand(K, 1), K); % 简化示意 centers = load_nodes(init_idx, :); for iter = 1:max_iter % 1. 根据当前站址做Voronoi分区 polys = bounded_voronoi(centers, bounds); % 2. 对每个分区计算负荷重心 new_centers = centers; for k = 1:K in_region = isinterior(polys(k), load_nodes(:, 1), load_nodes(:, 2)); nodes_k = load_nodes(in_region, :); loads_k = node_loads(in_region); if isempty(nodes_k) continue; end g = load_centroid(nodes_k, loads_k); % 3. 引力项:朝向重心移动 new_centers(k, :) = centers(k, :) + alpha * (g - centers(k, :)); end % 4. 斥力项:站与站之间保持安全距离 for i = 1:K for j = i+1:K dist_ij = norm(new_centers(i, :) - new_centers(j, :)); if dist_ij < safe_dist repel_dir = (new_centers(i, :) - new_centers(j, :)) / (dist_ij + eps); new_centers(i, :) = new_centers(i, :) + beta * repel_dir * (safe_dist - dist_ij); new_centers(j, :) = new_centers(j, :) - beta * repel_dir * (safe_dist - dist_ij); end end end % 5. 边界约束:站址必须在研究区域内 new_centers(:, 1) = min(max(new_centers(:, 1), region_x(1)), region_x(2)); new_centers(:, 2) = min(max(new_centers(:, 2), region_y(1)), region_y(2)); centers = new_centers; end这个循环跑完之后,我一般再叠加一个粒子群微调。原因很简单:围捕算法的特点是探索性强、空间分散性好,但局部精细搜索能力偏弱。粒子群在围捕得到的站址附近做精细搜索,能把目标函数再压低几个百分点。混合策略在多个随机场景下都比单独用一种算法效果好,具体数字后面实验部分说。
粒子群微调的代码比较常规,变量是站址坐标,目标函数是内层的Voronoi分区+潮流计算。粒子数取30,迭代30轮,惯性权重0.7,加速系数1.5。
4.4 潮流计算与目标函数评估
潮流计算是内层评估的核心环节。我手写了一个牛顿-拉夫逊法求解三相平衡配电网潮流,节点用极坐标形式。对300个节点来说,迭代在5次以内就能收敛,单次耗时在0.2秒左右。
关键实现点:
- 节点导纳矩阵 Y 用稀疏矩阵存储,避免全矩阵运算内存爆炸。
- 电动汽车充电负荷作为恒功率负荷(PQ节点)处理,但大规模接入时要修正为恒阻抗+恒功率混合模型,否则高渗透率场景潮流可能不收敛。
- 分布式光伏作为 PQ 节点(给定有功和功率因数),不是 PV 节点,因为配电网里光伏通常不参与调压。
潮流计算完成后,网损是直接可用的结果。投资费用部分我维护了一个station_config结构体,每个站的配置变量包括:变压器容量、充电桩数量、线路长度。这些和Voronoi分区的面积、负荷总量直接相关——面积决定了低压线路长度,负荷总量决定了变压器容量。
目标函数代码大致长这样:
function cost = evaluate_solution(centers, data) % 1. Voronoi分区 polys = bounded_voronoi(centers, data.bounds); % 2. 遍历分区,计算每个站的配置需求和费用 total_cost = 0; total_loss = 0; violation_penalty = 0; for k = 1:size(centers, 1) in_region = isinterior(polys(k), data.nodes(:, 1), data.nodes(:, 2)); nodes_k = data.nodes(in_region, :); loads_k = data.node_loads(in_region, :); % 站容量:分区峰值负荷 + 裕度 peak_k = max(loads_k) * 1.2; station_cost = station_investment(peak_k); % 投资等年值 % 线路费用:分区面积近似估算低压线路长度 area_k = area(polys(k)); line_cost = line_investment(sqrt(area_k) * 0.8); % 经验系数 total_cost = total_cost + station_cost + line_cost; end % 3. 潮流计算得到网损 [loss_mw, voltage_violation] = power_flow_with_violation(centers, data); total_loss = loss_mw * 8760 * data.elec_price / 1000; % 年网损费用(万元) % 4. 惩罚项 violation_penalty = voltage_violation * 1e6; cost = total_cost + total_loss + violation_penalty; end这里面有几个我踩过坑后确定的经验值:容量裕度系数取1.2,太保守会大幅抬高投资费用,太激进运行阶段容易过载;线路长度估算用sqrt(area_k) * 0.8这个经验系数,是在几个典型配电网络拓扑下拟合出来的平均值,做方案对比没问题,但如果要出严谨的工程结论,应该用实际路网数据做缓冲区分析。
5. 实验设计、参数调优与结果分析
5.1 对比实验怎么设计
为了验证“Voronoi+围捕”这套组合的有效性,我设计了四组对比:
- 方案A:只用K-means聚类站址,分区用Voronoi,不做围捕迭代。
- 方案B:Voronoi分区 + 粒子群搜索站址(没有围捕预搜索)。
- 方案C:Voronoi分区 + 围捕算法迭代 + 粒子群微调(本文完整方案)。
- 方案D:全区域单站方案,作为基准下限。
每种方案都跑10个随机种子,取平均值和标准差,避免单次随机结果带来误导。站点数量K分别取4、6、8三档,观察不同密度下的表现差异。
5.2 关键参数表与调参经验
先给一份我最终确定的参数表,供你复现时参考:
| 参数 | 取值 | 说明 |
|---|---|---|
| 研究区域范围 | 10km × 10km | 矩形区域 |
| 负荷节点数 | 300 | 随机生成,固定种子 |
| 电动汽车数量 | 5000 | 蒙特卡洛抽样 |
| 慢充功率 | 7kW | 居民区为主 |
| 快充功率 | 60kW | 商业区为主 |
| 引力系数 α | 0.35 | 调大收敛快,过大震荡 |
| 斥力系数 β | 0.15 | 调大站址更分散 |
| 站间安全距离 | 2.5km | 小于该值触发斥力 |
| 围捕迭代次数 | 30 | 实测20轮后基本收敛 |
| 粒子数/迭代数 | 30/30 | 粒子群微调配置 |
| 折现率 | 8% | 等年值折算 |
| 设备寿命 | 15年/20年 | 充电站/线路 |
调参过程中最值得说的一点是:α 和 β 的比例关系对结果影响非常大。α 太大,站址全被吸到负荷中心,围捕退化成K-means;β 太大,站址过于分散,远离负荷热点。我试过的组合里,α=0.35、β=0.15 是比较稳的搭配。还有一个细节:斥力计算要在引力更新之后做,如果先做斥力再做引力,站址会更容易震荡。
5.3 从仿真结果到规划结论
四组方案在K=6场景下的典型结果如下(年综合费用,单位万元):
| 方案 | 年综合费用均值 | 标准差 | 网损费用占比 |
|---|---|---|---|
| 方案A(K-means+Voronoi) | 2860 | 115 | 18.6% |
| 方案B(PSO+Voronoi) | 2750 | 85 | 17.2% |
| 方案C(围捕+PSO+Voronoi) | 2630 | 62 | 15.8% |
| 方案D(单站基准) | 3420 | — | 24.1% |
几个值得注意的发现:
第一,方案C比方案A年费用降低约8%,主要来自网损的下降和更合理的站点容量配置。围捕迭代后站址分布明显更均匀,各Voronoi区域的负荷规模标准差下降了约35%,这意味着没有出现“个别站超载、个别站闲置”的尴尬情况。
第二,方案C的标准差比方案B小,说明混合策略的稳定性更好。粒子群单独搜索站址时,搜索前期容易在某个局部区域徘徊,围捕预搜索相当于给粒子群提供了一个质量更高的初始种群,后面精细搜索的效率自然更高。
第三,K从4增加到8的时候,年综合费用并不是单调下降的——K=6比K=8低了约3%。原因很好理解:站多了投资费用上去,但网损下降幅度跟不上,边际效益为负。这类曲线做出来之后,决策者可以结合实际资金约束选择站点数量,而不是拍脑袋定。
5.4 可视化结果怎么看
我最后输出的图包括:
- 基础场景图:负荷点颜色按类型区分,电动汽车充电热力图用hexbin叠加。
- Voronoi分区图:每个分区用不同颜色填充,站址用星号标记,负荷重心用圆圈标记。围捕迭代过程中站址的移动轨迹用灰色细线画出,可以直观看到站址从初始位置往需求热点聚拢,同时彼此保持间距。
- 收敛曲线:横轴是迭代次数,纵轴是目标函数值,可以看到围捕阶段快速下降、粒子群阶段缓慢微调的两段式曲线。
- 方案对比柱状图。
Voronoi分区图是最值得仔细看的输出。我建议你输出每个分区的负荷总量和面积做成一个表格,这样能直观发现是否有分区面积很大但负荷很小的情况——如果有,说明该区域的站容可以调小,否则浪费投资。
另一个很重要的可视化是充电等待时间分布。我在每个分区内按充电需求到达率做了一个简单的M/G/k排队估算,画出全区域充电等待时间热力图。方案C下,最大等待时间比方案A减少了约20%,这个指标在向非技术背景的决策者汇报时比网损、潮流约束更有说服力。
6. 常见问题与排坑实录
6.1 Matlab实现中的高频报错
做一个问题速查表,这些都是我在复现过程中真实遇到过的:
| 问题现象 | 可能原因 | 解决办法 |
|---|---|---|
polyshape报错“顶层必须唯一” | Voronoi顶点共线或重合 | 调用simplify清理;站址加最小间距约束 |
| Voronoi区域出现NaN坐标 | 站址位于区域边界上 | 给站址加微小扰动或用unique去重 |
| 潮流迭代不收敛 | 电动汽车渗透率过高导致负荷模型失真 | 切换到恒阻抗+恒功率混合模型,或减小步长 |
| 优化迭代站址震荡剧烈 | α 取值过大 | 降到0.2-0.3,或加一个差分惯性项 |
| 粒子群重复搜索同一区域 | 初始种群多样性不足 | 改用Sobol低差异序列生成初始粒子 |
| 内存持续增长 | polyshape对象在循环里累积 | 每轮结束手动clear不需要的对象,或改用函数封装避免工作区污染 |
6.2 K-means初始化的坑与正确姿势
我前面提到围捕迭代前用K-means做初始化。这个选择背后有个非常现实的原因:Voronoi分区对初始站址非常敏感。如果你完全随机给初始站址,可能会出现初始时刻某个站点的区域内没有任何负荷点,它的重心没法算,围捕迭代就会异常。
K-means初始化之后,每个站点区域基本都有负荷了,围捕迭代就稳定了。如果你不想用K-means,也可以用K-means++的思想做最远点采样,保证初始站址尽可能分散开。实测下来,K-means初始化 + 围捕迭代的组合,在多数随机种子下都能得到质量稳定的解,波动系数在5%以内。
另外一个容易忽视的点是:不同K值要分别做初始化。K=4、K=6、K=8的结果不能直接从K=6的初始解插值出来,实测这样做的效果反而比重新初始化差,因为初始解空间的结构不同。
6.3 围捕迭代真正收敛了吗
这是我觉得最值得提醒的坑。只看目标函数曲线,你可能在第10轮就发现已经“看起来稳定了”,就提前终止迭代。但真的做下去,第15轮、第20轮的目标函数差异可能只有0.5%,但对最终站址布局的影响却很明显。
我的建议是,不要只盯着目标函数值,要看分区后的负荷均衡度和站间距离。这两个指标的收敛速度比目标函数慢得多。目标函数稳定不代表分区结构稳定,而分区结构才是规划方案真正需要确定的东西。
我做了一个检查函数,每轮迭代后输出所有分区的负荷标准差,直到标准差变化率小于1%才判定收敛。视觉上,你会看到站址在第21轮时还在缓慢移动,虽然目标函数已经非常平稳。这个“先看负荷分布、再看目标函数”的收敛判断准则,我强烈建议你也用上。
6.4 关于仿真时长和性能优化
300个节点、6个站点、围捕30轮 + 粒子群30轮,每一轮都要做Voronoi分区、负荷分配、潮流计算,运行时间大概是15分钟。如果你在一个大区域、1000个节点、10个站的场景下跑,时长直接翻几倍。
我做了三件提速的事:
- 用
parfor并行评估粒子群里的每个粒子。粒子间完全独立,天然适合并行。6核机器实测加速约4.5倍。 - 站址更新只对围捕迭代中的分区域做面积和负荷重心的增量更新,不对全局重新计算。这个优化在线性代数层面减少了很多重复计算。
- 把潮流计算的稀疏矩阵因子化结果缓存起来,只有节点负荷变化超过阈值时才重新分解矩阵。因为规划阶段的潮流变更主要是负荷水平整体变化,拓扑基本不变。
跑完一遍之后我还做了一次敏感性分析:把电动汽车渗透率从10%逐步提升到50%,观察Voronoi分区的边界如何移动、站址重心如何向西侧偏移。这类分析输出的结论比单一方案的数值结果更有参考价值,因为规划者关心的是“未来不同发展情景下方案怎么适应变化”。
7. 延伸思考:这套方法还可以怎么用
这个项目做下来,我对“Voronoi图+围捕算法”这个组合的理解深了不少。其实这套东西完全可以迁移到类似的时空优化问题上,不一定非要限定在配电网:
- 共享单车调度区域划分:把调度站当站址、骑行需求热点当负荷,用同样的双层框架去优化调度站布局。
- 物流快递网点规划:网点选址和配送范围划分,本质上和充电站选址是同一个数学结构。
- 分布式储能规划:储能位置和容量的配置跟充电站同理,只是目标函数里的收益模型不同。
- 多功能车巡游路线规划:巡逻点布局、无人机机巢部署,只要有“设施选址+服务范围划设”的问题,这套方法都值得一试。
迁移的时候只需要替换两件事:一是负荷/收益模型,二是目标函数里的成本和约束定义。辛算法框架、Voronoi分区、围捕迭代和粒子群微调的过程完全复用。
最后再分享一个我个人的体会:这类跨学科组合的项目,难点从来不在单个算法本身——Voronoi图、粒子群、潮流计算,任何一个单独拿出来都有大量成熟资料。真正花时间的是把它们拼起来之后,接口处的细节处理。比如Voronoi裁剪边界的数值稳定性、站址更新和分区之间的顺序依赖、优化算法与评估模型之间的收敛判定。这些接口问题在论文里通常一句话带过,落在代码里却是实打实要一个一个解决的坑。
我第一次跑通这个项目的时候,卡在Voronoi边界裁剪上将近一周,最终搞定的那一刻反而没有太多兴奋感——因为后面等待着的是参数调优、对比实验、结果分析这些更琐碎但也更重要的事情。如果你也正在复现类似的算法组合,遇到边界处理、收敛震荡、结果不稳定这些问题,希望这篇文章能帮你少走一些弯路。如果后续你跑出了和我不同的参数结论,欢迎交流差异背后的原因,那通常意味着你的场景约束里藏着新的有价值的信息。