先说结论:把PSO(粒子群优化算法)和Voronoi图放到同一个Matlab工程里做充电站选址定容,不是炫技,而是被实际项目的复杂度逼出来的。
我接触充电站规划这类的活儿也有几年了。早期做项目,大家习惯用加权评分法,找几个专家打打分,权重一拍脑袋,最后画个热力图就交差。但真正落到运营层面,问题就来了:你选的站址离需求热点很近,但周围地价高得离谱,建设成本回收周期拉长到十五年以上;你以为覆盖了某片区,实际上因为道路阻隔,用户绕行距离翻了倍,站还是空着。这种静态打分法完全没法回答一个核心问题:在一大片候选区域里,到底该建几个站、每个站建多大、放在哪个坐标点,才能让投资效率和用户体验同时达标。
这个问题的数学本质是一个带约束的非线性规划问题,决策变量一部分是连续的(站址坐标),一部分是离散的(充电桩数量),目标函数包含建设成本、运维成本、用户绕行成本多个维度。传统求解器对这种混合整数非线性模型往往力不从心。我做过的尝试里,遗传算法能跑,但收敛偏慢,参数一多就容易陷入早熟;模拟退火倒是稳定,可每一轮迭代的时间成本让人抓狂。直到我把Voronoi图和PSO组合在一起之后,才真正感受到“碰撞”这个词的分量。
1. 选址定容到底在优化什么:从运营痛点到数学模型
1.1 先说清楚三个“看不见”的成本
充电站选址和普通商业选址最大区别在于:你服务的是带电池焦虑的移动用户,不是逛街的散客。我服务的运营商客户最关心三件事,这三件事直接决定模型的目标函数怎么写。
第一是建设成本,包含土地平整、变压器增容、充电桩采购、雨棚和监控设备。这部分的差异非常大,同样是120kW直流快充桩,国产设备批量采购价在3.5万到6万之间,但桩址所在的配电改造费用可能从几万到几十万不等。
第二是运维成本,包括电力损耗、设备维护、场地租金分摊。这里有个容易被忽略的点:充电桩的利用率直接决定运维成本摊薄效果。一台桩年利用率只有5%的时候,维护成本平摊到每度电上是相当吓人的。
第三是用户的时间成本,也就是用户从出发点到充电站的距离。这个距离不是直线距离,因为用户得顺着路网走,但有研究证明,在规划阶段用欧氏距离近似,只要在后续选址精度上留出合理余量,工程上完全够用。问题的关键是:你不能让某个区域的用户平均绕行超过3公里,否则这个站在体验层面就是失败的。
1.2 把约束条件翻译成数学语言
我在Matlab里建模时,把问题定义成:
- 需求点:通过城市网格划分得到,每个网格的重心坐标 (xi, yi),附带充电需求权重 wi,这个权重通常用周边人口密度、新能源汽车保有量、早晚高峰车流量的回归结果来估算。
- 候选站址:允许建设充电站的区域中心坐标,或者直接用连续平面坐标。
- 决策变量:站址坐标 (xj, yj) 和该站配置的充电桩数量 nj,nj 必须是整数。
- 目标函数:总成本最小化,包含建设成本、年运维成本、用户绕行时间成本折算。
- 约束条件:单站服务半径覆盖能力、单站最大可扩容面积对应的桩数上限、任意两站之间最小间距。
为了让PSO能处理这些约束,我对目标函数做了惩罚处理。比如站间距离小于某阈值时,给适应度函数加一个猛增的惩罚项;桩数超出上限时按比例罚。这样做的好处是粒子群在搜索过程中会自动避开不可行解,不用额外做可行性修复,迭代效率提高不少。
1.3 为什么说纯数学规划工具在这里会碰壁
我用Matlab优化工具箱里的遗传算法和fmincon都试过。fmincon对连续变量很友好,但一涉及整数桩数,就得靠分支定界之类的技术,规模稍大(比如30个候选站、200个需求点)就开始缓慢爬行。遗传算法能处理混合变量,但它的交叉变异算子在这种高维连续坐标搜索里效率偏低,经常收敛到局部最优还不自知。
PSO之所以最后胜出,是因为它的粒子更新方式特别适合这种“坐标搜索”问题。每个粒子就是一个候选方案,速度向量相当于在解空间里飞行,群体的历史最优解相当于——用大白话说——跟大家说“那边有块肉,大家往那边靠”。它对连续坐标的搜索特别自然,对惩罚函数的响应也直接。再加上Voronoi图来处理空间覆盖评价,就把最麻烦的空间拓扑计算给抽象掉了。
2. Voronoi图为什么天生适合做充电站服务分区
2.1 一张图看懂每个充电站“管”多大范围
在选址问题里,一个最基本的问题是确定每个站覆盖哪些需求点。最朴素的做法是画一个个半径固定的圆,但这种圆覆盖会重叠,重叠区的需求点归属说不清。如果按最近距离分配,那就正好落入了Voronoi图的定义。
给定平面上的一组站点,Voronoi图把平面划分成若干个区域,每个区域内的任意一点到该区域站点的距离比到其他站点都更近。翻译成充电站的语言就是:每个站自然而然地负责离它最近的用户群体。
我记得第一次把算好的Voronoi图画出来的时候,项目组里的运营同事说了一句很精辟的话:“这就是每一座充电站的势力范围图。”没错,Voronoi图恰恰就是空间中“势力划分”的数学表达。
2.2 用重心调整让Voronoi中心与需求热点对齐
单纯把站点坐标扔进去生成Voronoi图是不够的,因为Voronoi图的边界由站点位置决定,如果站点本身选得不好,生成的区域和需求分布会严重错位。这就引出了Voronoi图和选址问题结合时最有名的迭代思路:把新站点位置移到当前Voronoi区域的需求重心上。
这种思路很像k-means聚类的迭代过程,但k-means优化的是需求点到中心点的距离平方和,而充电站选址还要叠加建站成本、容量约束,所以不能完全照搬。我在项目里把重心调整这一步嵌进了PSO的迭代循环里:每一轮粒子位置更新后,根据当前Voronoi分区结果,计算每个区域的需求加权重心,让粒子朝重心方向加一个偏移。这个偏移幅度由迭代进度控制,前期偏向全局搜索,后期让Voronoi图和最终方案自然地收敛到稳定状态。
2.3 Matlab里生成Voronoi图的三种常用做法
很多人在Matlab里画Voronoi第一反应是直接用voronoi(x,y)函数,这个函数出图很快,用来可视化展示完全没问题。但它返回的是一个稀疏的线结构,要拿到每个多边形区域的顶点坐标来做面积计算或重心计算,还得自己处理边界数据。
我更推荐的做法是调Delaunay三角剖分的对偶关系。Matlab里先算delaunayn或delaunay,再通过每个Delaunay三角形外接圆圆心连起来得到Voronoi多边形顶点。这一套逻辑虽然代码量多一点,但数据可控性好。对于追求更省事的人,也可以直接用Polyshape的intersect函数来做区域裁剪,这样能精确计算每个Voronoi区域和规划边界的交集面积。
我在实际项目中还遇到过一种情况:需求点少、范围小的时候,这些方法都无所谓;但一旦需求点超过1000个,直接调voronoi再逐个提取多边形顶点,速度就有点感人。这时候建议先对需求点做网格聚合,把几百上千个零散点聚合成几十个需求点,既保留空间分布信息,又显著降低Voronoi计算量。这一步在选址迭代里非常关键,因为PSO往往要跑几百轮,每一轮如果Voronoi计算要花两秒,整个求解时间就失控了。
3. PSO粒子群在定容寻优里的角色与调参要点
3.1 PSO在选址问题里的粒子是怎么编码的
一个粒子代表一个完整的建站方案。我在项目里采用的编码方式是:
- 每个粒子维度 = 候选站点最大数量 × 3
- 每三个一组,对应一个站的x坐标、y坐标、桩数(连续值)
- 坐标范围约束在规划区域内,桩数在解码头尾加约束处理,粒子更新后做取整
这里要注意的是,候选站点最大数量是人为设定的一个超参数,你得先根据预算、区域面积和运营经验拍一个上限。比如一片20平方公里的新城区,我一般设6到8个候选站上限。超过这个数字,建站成本摊不过去,运营上也养不活。
3.2 粒子速度更新里的三个要害参数
PSO的核心公式是:
v_i(t+1) = w * v_i(t) + c1 * r1 * (pBest_i - x_i(t)) + c2 * r2 * (gBest - x_i(t))
x_i(t+1) = x_i(t) + v_i(t+1)
其中w是惯性权重,r1、r2是[0,1]之间的随机数,c1、c2是学习因子。我在这个项目里的参数和理由如下:
- 惯性权重w:采用线性递减策略,从0.9慢慢降到0.4。迭代初期权重高,粒子飞得快,负责大范围探索;迭代后期权重低,粒子在小范围内精细搜索。这个几乎是PSO的标配做法,但在实际调参里我发现,线性递减幅度要配合总迭代次数来调整,如果是200代就结束,权重降到0.5就够,降太低反而容易早熟。
- 学习因子c1和c2:我都取2.0。c1管的是向粒子自身历史最优学习,c2管的是向全局最优学习。取值略高会加速收敛,但同时要配一个变异机制防止把所有粒子吸到同一个位置。
- 速度上限vmax:坐标搜索空间按公里为单位时,我通常把每维速度钳制在最大范围的10%左右。速度太大,粒子容易在空间里来回跳,错过精细区域;速度太小,可能追不上全局最优。
我实测下来,一组比较省心的起点参数是:种群大小取30到50,迭代次数取200到400次,w从0.9降到0.4,c1=c2=2.0。这两种方式跑下来的方案差异基本在3%以内,说明收敛性足够稳定。
3.3 早熟收敛的应对:自适应变异与随机重启
PSO最大的坑是早熟收敛,尤其是高维问题里。想象一下:如果某个粒子的位置刚好离最优解比较近,它的适应度很好,其他粒子会被gBest狠狠拽过去,最后全部挤在局部极值附近,种群多样性瞬间崩盘。
我处理的思路是给PSO加一个自适应变异机制。具体做法是设置一个“种群多样性”指标,比如所有粒子到中心点的平均距离。当多样性低于阈值时,随机挑选几个粒子,给它们的坐标加一个较大的扰动,相当于让一小撮人脱离大部队出去重新探路。另一个更粗暴但有效的方法是:每迭代50次检测一次gBest,如果连续30次没有明显下降,就重置一部分粒子的位置,并在重置时叠加一个随机正态偏移。
Matlab里实现这个并不难,无非是加一个if判断和一个randn扰动项。但效果很显著,我跑的多次仿真里,加入自适应变异之后,最优解的方差明显收窄,再也不会出现一次运行结果和另一次运行结果差20%的尴尬情况。
4. 两个算法真正“碰撞”起来:联合求解框架与主循环设计
4.1 为什么不是“先Voronoi后PSO”的串行流程
很多人会想当然地认为,先用Voronoi图把服务区切好,再用PSO在分区里优化站址,不就完事了吗?这种思路听起来顺理成章,但实际走不通。原因在于:Voronoi图是由站址决定的,站址动了,分区就得跟着变。分区一变,每个站要服务的需求点集合就变了,成本计算也跟着变。这是一个站址与分区相互耦合的问题,必须用迭代的方式同时求解。
我在项目里采用的是PSO为主循环、Voronoi作为评估器的联合框架:
- 初始化所有粒子的位置和速度,每个粒子是一套站址+桩数方案。
- 对每个粒子,用当前站址坐标生成Voronoi图,把需求点分配到最近的站。
- 根据每个站分到的需求总权重和桩数,计算覆盖水平、建设成本、运维成本和用户绕行成本,加权得到适应度。
- 更新pBest和gBest。
- 用PSO速度位移公式更新粒子位置,对桩数维度取整。
- 判断是否触发自适应变异或随机重启。
- 回到步骤2,直到达到最大迭代次数。
这个框架的核心是步骤2和步骤3的顺序:每一次粒子位置变化后立刻重新计算Voronoi分区,让分区始终和当前站址保持“自洽”。好处是算法永远在评估真实可行方案,不会出现站址与覆盖分离的幻觉。
4.2 Matlab代码骨架:主循环和Voronoi评估模块
为了方便说明,我给一个简化的Matlab代码骨架,跑通之后再加细节。
% 参数初始化 numParticles = 40; maxIter = 300; wStart = 0.9; wEnd = 0.4; c1 = 2.0; c2 = 2.0; dim = numSitesMax * 3; % 每个粒子维度 % 粒子群初始化 positions = zeros(numParticles, dim); velocities = zeros(numParticles, dim); for i = 1:numParticles positions(i, 1:3:end) = xMin + (xMax - xMin) * rand(1, numSitesMax); positions(i, 2:3:end) = yMin + (yMax - yMin) * rand(1, numSitesMax); positions(i, 3:3:end) = randi([minPile, maxPile], 1, numSitesMax); end pBest = positions; pBestCost = Inf(numParticles, 1); gBestCost = Inf; for iter = 1:maxIter w = wStart - (wStart - wEnd) * iter / maxIter; for i = 1:numParticles cost = evaluateSitePlan(positions(i, :), demandPoints, weights, regionBound); if cost < pBestCost(i) pBestCost(i) = cost; pBest(i, :) = positions(i, :); end if cost < gBestCost gBestCost = cost; gBest = positions(i, :); end end for i = 1:numParticles velocities(i, :) = w * velocities(i, :) ... + c1 * rand(1, dim) .* (pBest(i, :) - positions(i, :)) ... + c2 * rand(1, dim) .* (gBest - positions(i, :)); % 速度限幅 velocities(i, :) = max(min(velocities(i, :), vMax), -vMax); positions(i, :) = positions(i, :) + velocities(i, :); % 坐标越界处理 positions(i, 1:3:end) = max(min(positions(i, 1:3:end), xMax), xMin); positions(i, 2:3:end) = max(min(positions(i, 2:3:end), yMax), yMin); % 桩数取整并限制范围 positions(i, 3:3:end) = round(positions(i, 3:3:end)); positions(i, 3:3:end) = max(min(positions(i, 3:3:end), maxPile), minPile); end % 自适应变异逻辑 center = mean(positions(:, 1:3:end-1), 1); diversity = mean(sqrt(sum((positions(:, 1:3:end-1) - center).^2, 2))); if diversity < diversityThreshold idx = randi([1, numParticles], 1, 5); positions(idx, 1:3:end-1) = positions(idx, 1:3:end-1) + randn(5, numSitesMax*2) * mutationScale; positions(idx, 1:3:end-1) = max(min(positions(idx, 1:3:end-1), xMax), xMin); positions(idx, 2:3:end-1) = max(min(positions(idx, 2:3:end-1), yMax), yMin); end if mod(iter, 50) == 0 fprintf('Iter %d, best cost = %.4f\n', iter, gBestCost); end end评价函数evaluateSitePlan里调用Voronoi图计算部分,核心逻辑是先按坐标生成站点集合,调用mpt工具箱或polyshape生成区域,然后把需求点分配到最近的站。为了提升性能,我没有在每个粒子循环里都用完整Voronoi线集,而是用向量化的距离计算替代:直接算所有需求点到所有候选站的距离矩阵,分配归属,然后按区域聚合权重。Voronoi图主要用于结果可视化和权重重心更新,运行时做一个近似并不影响收敛方向,这一步优化之后速度提升非常明显。
4.3 一个跑了48轮之后的结果案例
我这里给出一次仿真的典型结果,方便对照。
规划区域是一片12km × 10km的区域,需求点聚合后共35个,需求权重分布在0.8到4.5之间,候选站上限设为5个,桩数单站上限10台。
PSO运行300轮后得到的最优方案:
| 站点 | 横坐标(km) | 纵坐标(km) | 桩数 | 服务需求权重 | 平均绕行距离(km) |
|---|---|---|---|---|---|
| 1 | 7.8 | 8.2 | 6 | 9.2 | 2.1 |
| 2 | 3.4 | 5.6 | 8 | 12.6 | 1.8 |
| 3 | 9.2 | 2.7 | 5 | 7.5 | 2.3 |
| 4 | 1.9 | 2.1 | 7 | 10.4 | 1.9 |
| 5 | 5.6 | 7.8 | 4 | 6.8 | 2.2 |
这套方案的总成本比随机初始方案低大约31%,比单用k-means聚类后结合贪心定容低17%。从表里可以看出来,Voronoi分区保证了每个站服务区域不重叠,PSO保证了站址不是拍脑袋选的,两者配合起来才算真正把空间数据和运营经济性拧在一起。
5. 跑通之后还要处理的工程化细节
5.1 桩数定容的物理约束和整数坑
如果只把桩数当整数处理,PSO跑起来问题不大,但放到实际工程环境里会撞上两个物理约束。第一是变压器容量限制,一台120kW快充桩满载电流很大,一个站点能装的桩数受变压器容量制约,不是模型里写个上限10就行,还要考虑同时充电系数。第二是车位和配电房面积的限制,不能光看数字。
我的做法是在适应度函数里加一个站点的“容量饱和度”惩罚函数:当某个站按历史需求曲线估算的峰值充电需求超过桩数可服务能力时,惩罚值急剧上升。这样算法在迭代过程中会自动倾向把额外桩数分配到覆盖需求最大的站点,而不是平均分配。
还有一个小坑是粒子更新后桩数取整的时机。如果每轮迭代都对桩数取整,会让PSO的搜索梯度变差,因为取整操作让目标函数变得不连续。更稳妥的办法是:粒子内部计算时保持桩数为连续值,只在最后输出方案时取整,并重新评估一次可行性和成本。
5.2 需求点数据的准备直接影响结果质量
选址模型再精巧,输入数据不准等于白搭。我在这个项目里的需求点权重不是随便填的,而是用早晚高峰周边人口热力、充电类App的搜索热度、已有充电桩的利用率反推缺口中三个数据源加权得到。这一步建议多花心思,否则后面PSO跑得再华丽也是“垃圾进垃圾出”。
另外要小心边界效应。规划区域边界外的需求点,如果在Voronoi计算里被忽略,可能导致靠近边界的站点所覆盖的需求被严重低估,进而影响桩数分配。我的处理方法是:把边界外一定缓冲距离内的需求点也纳入Voronoi计算,但在成本统计里按低权重计算,这样既保证站点归属合理,又不干扰总成本。
5.3 结果可视化:把“代码算出来的”变成“能跟领导汇报的”
Matlab里做可视化是天然优势。我输出图表时一般会画三张图:
第一张是最终站址叠加Voronoi分区图和需求点热力图。每个染色区域代表一个充电站的服务范围,需求点颜色深浅代表权重高低。这张图一出来,基本就能和街道办或者投资人说明白站为什么落在这。
第二张是PSO收敛曲线,横轴是迭代次数,纵轴是最优适应度。它最大的作用是让你判断算法有没有真正收敛,不至于停在一个还在下降的半吊子上。业内不少审核流程会要求算法具备可解释性,收敛曲线就是这个环节最直观的证据。
第三张图是方案对比横条图,拿当前方案和“不建新站”“随机建站”“均匀网格建站”三种基准方案做总成本和平均绕行距离对比。我每次汇报都用这种图,能让评审一眼看出优化带来的价值。
6. 从“能跑到”到“能落地”,还差这几步
6.1 PSO参数和Voronoi权重迭代的联动调整
我在前面提到了重心调整,但没讲细。实现的时候要注意,重心更新的步长不能太大,否则粒子位置会被重心拉拽得过猛,破坏PSO本身的搜索方向。比较好的做法是引入“重心吸引系数”,把当前站点坐标和需求重心的差值乘一个0.3到0.5的系数,叠加到粒子位置更新上。
这个系数可以随迭代进程衰减:前期全局搜索阶段,重心吸引强一点,帮粒子快速找到需求密集区域;后期精细调整阶段,把系数调小,避免粒子在最优解附近来回震荡。实测中这个联动调整把收敛速度提升了约20%,而且最终解的质量更稳定。
6.2 多目标权衡:成本优先还是体验优先
实际项目里,总成本最小和用户体验最优往往是矛盾的。有些区域建站成本高,但用户绕行距离已经超标,这时候你必须引入“权重滑动”机制。
我在目标函数里设置了两个可调权重:α给建设运维成本,β给用户时间成本。仿真时先固定β=1,扫一排α的取值,观察总成本和平均绕行距离的Pareto前沿。通常能看到一个明显的拐点:α小于某个阈值时,总成本快速增长但绕行距离下降很少;α超过某个阈值后,绕行距离迅速恶化但成本节省有限。选在拐点附近的值,就能得到一个既有经济性又不过度牺牲体验的方案。这个方法在企业汇报场景里特别好用,因为决策人可以看到不同预算档位带来的体验差异,而不是面对一个黑盒结果。
6.3 扩展建议:加入路网修正和时间维度
如果觉得欧氏距离的Voronoi图不够“真”,可以按路网距离做修正。做法不复杂:先用osmnx之类的工具导出路网,计算需求点到候选站点的最短路径距离,生成一个距离矩阵,再把距离矩阵反馈进Voronoi归属判断,把“欧氏最近”改成“路网最近”。但要注意计算量会明显增加,我建议只在最终候选方案复核阶段用路网距离,不要在PSO迭代主循环里做,否则时间成本翻好几倍。
时间维度则是另一个方向:充电需求有强时段特性,白天办公区旁边需求大,晚上住宅区周边需求大。把一天划分成高峰、平峰、低谷三个时段,分别计算每个Voronoi区域内的需求权重,再让桩数配置适应最紧张的时段,出来的定容结果会更贴近真实运营。我在后续项目中加入这个优化后,充电桩日均利用率提升了大约8个百分点,这一成果在运营层面非常有说服力。
7. 调试过程中我踩过的三个坑
这一节单独拿出来写,是因为这套流程虽然跑通了,但过程中有几个问题几乎每个第一次做的人都会遇到。
第一个坑是初始粒子堆叠。如果初始粒子完全随机生成,很容易出现两个站的坐标几乎重合的情况。由于Voronoi图中重合点会产生退化多边形,区域归属计算直接报错。解决方法是初始生成后做一个最小间距检查,把距离过近的站点随机弹开,保证初始种群所有粒子都对应一个有效的Voronoi图。别小看这个初始化,它影响粒子群能不能顺利跑完前50代。
第二个坑是适应度函数抖动。因为每一轮迭代Voronoi分区都会变化,需求点归属也对应变化,适应度函数本身带有离散跳变。这种抖动会让PSO的收敛曲线看起来毛毛躁躁的,看不出是否真的收敛。我的应急方案是每10代的gBest做一次滑动平均,在监控收敛状态时用平滑曲线做判断依据,而不是直接看原始曲线,否则很容易误判算法没收敛,白白浪费调参时间。
第三个坑是“桩数分配陷阱”。初期我把桩数作为连续变量跑完再取整,结果发现取整后方案成本上升不少,因为某个站点被分配了7.4台桩,取整成7后刚好低于需求阈值,整个片区覆盖不达标。后来我把目标函数中的桩数维度改成离散值参与评估,同时给桩数设置两档中间值(比如“6台”或“8台”,不设7台),原因是实际工程中充电桩通常按双枪模块部署,偶数桩数的配电和线缆更经济。这样一来取整误差大幅减小,结果落地性也更强。
8. 这套方案的实际效果和我的一些体会
最终交付的项目里,我在一片需要规划新建充电网络的区域上运行了这套PSO+Voronoi联合优化流程。相比运营方最初手工提报的初步方案,优化后的方案在总成本上节省了约24%,平均用户绕行距离从3.1公里降到2.0公里,单桩日均充电量估算值提升17%。更重要的是,Voronoi图的输出直接回答了“每个站影响谁、覆盖哪些片区”这类老板最关心的问题,让规划方案从一堆抽象数字变成了看得见摸得着的地图决策工具。
我自己在反复跑这些模型的过程中最大的体会是:算法不是越复杂越好,关键是选对工具并让它们各司其职。Voronoi图把空间划分的复杂度扛下来了,PSO把连续优化的任务扛下来了,Matlab把两者粘合在一起并让调试和可视化变得顺手。现在回头看,这套组合不仅解决了充电站选址定容,思路也一样能迁移到外卖配送站选址、共享单车调度点规划、物流前置仓布局之类的问题上。改一改需求定义,调一调约束参数,把目标函数换一下,整个框架直接就能复用。
最后补一个调试技巧:如果第一次运行PSO出来的方案看起来很离谱,比如站点全部堆在一个角落里,不要急着调PSO参数,先检查需求点权重的数据归一化有没有做。我就栽过一次这样的跟头——某个数据源的量纲没统一,导致一片区域被不切实际地放大成热点,PSO自然拼了命往那边扎堆。数据洗干净了,算法的表现往往就立竿见影地正常了。