简介面向电动汽车充电站规划场景的Matlab程序包以规划期充电站总成本投资、运行维护与网损费用之和最小为目标在相关约束下构建电动汽车充电站最优选址定容数学模型从34个候选位置中优化选取7个站址。包内共4个Matlab脚本压缩包仅5KB包含主程序、基于Voronoi图的服务区域划分、成本计算等核心模块代码注释清晰便于初学者快速掌握建模思路。目前已有614人学习下载适合电力系统、交通电气化方向的学生用于选址定容问题的算法验证与课程设计。通过运行程序可完整复现候选点筛选、Voronoi划分、总成本核算的优化流程有助于理解充电站规划中的空间分析与数学规划方法也可作为相关论文或竞赛作品的Matlab参考实现。1. 为什么充电站选址要当成「34 选 7」的组合优化问题电动汽车充电站选址定容这个 matlab 程序要处理的是规划阶段最头疼的问题备选名单里躺着 34 个候选站址但预算和电网接入条件摆在那里只能选出 7 个来建还得顺带把每个站的容量定下来。凭感觉挑人流密集处建站往往忽略了配电网损和后期运维等到运营时才发现成本失控。这个程序把选址定容建模成组合优化问题用 Voronoi 图划分每个站的服务区目标函数是规划期内总成本投资运行维护网损最小。跑完能直接看到哪 7 个点入选、容量如何分配。代码由 main.m、VoronoiT.m、VorCostCDEV.m、VoronoiArea.m 组成注释清晰适合刚接触充电站规划的初学者和需要做课设的读者。2. 数学模型与成本构成总成本 网损怎么进目标函数2.1 目标函数拆解摘要里的描述已经把优化目标说清楚了规划期内充电站的总成本投资、运行和维护成本和网损费用之和最小。这里的关键是“之和”因为一次性投资和年度网损费用单位不同不能直接相加。常见做法是把投资成本按设备寿命折算成等年值运行维护成本按投资额的比例估算网损费用则基于典型日负荷曲线乘电价算出年费用这样三项才能放在同一个目标函数里。% 目标函数示意三项成本加总 function C totalCost(x, data) idx find(x 1); % 被选中的站点编号 Cinv sum(data.InvCost(idx)); % 投资成本已折算为等年值 Com sum(data.OmCost(idx)); % 年运行维护成本 Closs sum(data.Loss(idx)); % 年网损费用 C Cinv Com Closs; end这段代码把决策变量 x34 维 0-1 向量映射到总成本。data.InvCost、data.OmCost、data.Loss 是 34 维列向量分别代表每个候选点若建站时对应的年化投资、年运维费用和网损费用。注意这里只是示意结构实际程序里这些值是在 VorCostCDEV.m 中根据站址和所属负荷动态算出来的而不是预先给定的常数这样做的原因是网损费用与 Voronoi 分区结果强相关站选得越偏、负荷距离越远网损就越高。2.2 约束条件34 选 7 不是简单数数去掉约束的组合优化没有实际意义充电站选址至少要满足四类约束。第一是数量约束恰好选 7 个站即 sum(x) 7。第二是容量约束每个站点的装机容量要在合理区间 [Pmin, Pmax] 内同时该站 Voronoi 分区内所有负荷需求之和不小于容量下限、不超过容量上限否则要么利用率太低要么高峰期排队到马路对面。第三是服务半径约束负荷点到最近站的距离不能超过 Rmax超过的话用户会流失这个约束在分区后统计最大距离来校核。第四是站间最小距离约束防止两个站建在同一个路口互相抢客。% 约束条件伪代码 sum(x) 7 % 恰好 7 个站点 Pmin Pcap(idx) Pmax % 容量上下限 maxDist(voronoiCell(i)) Rmax % 服务半径 dist(site(i), site(j)) Dmin % 站间最小距离这些约束在 matlab 程序里很少全部写成等式约束大多数时候是加到目标函数里做惩罚项原因后面 4.2 会讲。Pcap 是每站的容量变量它既受到 Pmin/Pmax 限制也反过来决定 VorCostCDEV.m 里投资成本的大小所以容量和选址实际上是耦合在一起的这也是“定容”和“选址”不能拆开算的原因。2.3 变量与参数表符号含义说明N候选点数量本例 N34K计划建站数本例 K7x34 维 0-1 决策变量1 表示在对应候选点建站C_inv投资成本含土地、设备、配电设施折成等年值C_om运行维护成本常按投资额的一定比例估算C_loss网损费用与负荷空间分布、到站距离相关Pcap站点容量同时影响投资成本和可服务负荷量Rmax最大服务半径由用户接受度和项目要求决定表格里的参数在 main.m 里通过 data 结构体传入实际取值来自输入数据文件。需要说明的是不同文献对等年值系数的取法差异很大有的用 8% 折现率、10 年寿命有的用 12 年这会影响最终选站结果。我一般会先按题目给定值跑通再做一个敏感性分析看折现率变化是否会导致最优方案跳变。3. Voronoi 图分区与 Matlab 核心文件拆解3.1 VoronoiT.m站与站之间的势力范围Voronoi 图又叫做泰森多边形给定一组站址后平面会被划分成若干多边形每个多边形内的任意点到对应站址的距离都小于到其他任何站址的距离。对充电站来说这恰好模拟了用户“谁近就去谁那充电”的行为因此某个多边形内的负荷就应该由该站承担。VoronoiT.m 这个文件在程序里承担的就是这个划分工作它把站址坐标转换为每个站的多边形顶点。matlab 内置了 voronoin 函数可以返回顶点和单元索引VoronoiT.m 大概率是对它的一层封装同时处理边界裁剪和 Inf 顶点。常见写法如下% VoronoiT.m 的常见封装形式 function [cellVerts, cellIdx] VoronoiT(siteXY) [V, C] voronoin(siteXY); % V 为顶点坐标C 为单元顶点索引 nSite size(C, 1); cellVerts cell(nSite, 1); for i 1:nSite cellVerts{i} V(C{i}, :); % 取出第 i 个站的多边形顶点 end cellIdx C; end参数说明siteXY 是 34 行 2 列的坐标矩阵V 是所有 Voronoi 顶点的坐标C 是一个 cell 数组C{i} 是第 i 个站的顶点索引序列。需要特别注意的是位于凸包边界上的站其多边形会延伸到无穷远V 中出现 Inf导致后续 polyarea 计算得到 NaN。我一般会在 VoronoiT 之后加一步裁剪把多边形限制在规划区域范围内否则整个目标函数都会变成 NaN优化直接崩掉。3.2 VorCostCDEV.m单站成本怎么算这个文件名拆开看是 Voronoi、Cost、CDEV 的组合CDEV 大概指充电设备或电动汽车充电设施。它的输入应该是站址坐标、站点容量、该站所辖负荷点信息输出是该站的总成本包括投资、运维和网损。% VorCostCDEV.m 计算单个充电站的年成本 function cost VorCostCDEV(x, y, Pcap, loads, price, life) dist sqrt((loads(:,1) - x).^2 (loads(:,2) - y).^2); invCost 0.08 * (500 1200 * Pcap); % 年化投资估算 omCost 0.03 * invCost; % 运行维护按 3% lossCost sum(dist .* loads(:,3)) * price * 8760 / 1000; cost invCost omCost lossCost; end这里用负荷矩功率乘距离来近似网损费用是选址阶段常用的一种简化手段不需要建立完整潮流模型。loads 的第三列是各负荷点的充电功率price 是单位电价8760 是一年小时数除以 1000 是单位换算。实际项目里如果规划区域电网结构已知可以用 DistFlow 潮流计算替换这一行精度会更高但计算时间也会随之增加。这个函数的返回值就是 2.1 节目标函数中单个站的 C_inv C_om C_lossmain.m 需要把所有选中站的结果加起来。3.3 VoronoiArea.m用多边形面积校核容量VoronoiArea.m 解决的问题是拿到 VoronoiT 输出的多边形后怎么得到面积并据此校核容量。如果各负荷点的功率密度已知面积乘密度就是该站服务区内的总负荷如果负荷点离散则需要判断负荷点是否落在多边形内再累加功率。% 计算单个 Voronoi 单元的面积 function area VoronoiArea(cellVerts) if any(isinf(cellVerts(:))) area NaN; % 未裁剪时容易触发 else area polyarea(cellVerts(:,1), cellVerts(:,2)); end endpolyarea 是 matlab 内置函数直接输入顶点坐标就能返回多边形面积。上一层中如果 cellVerts 含 Inf这里返回 NaN进而导致 VorCostCDEV 里的负荷统计失效所以在使用前一定要先做边界裁剪。我见过不少初学者在 matlab 论坛上问“为什么画出 Voronoi 图有射线”其实就是没处理边界这不会让程序报错但会让面积计算静默出错。3.4 main.m 主流程主程序的任务是把上面几个函数串起来。读取数据后初始化一组满足 sum(x)7 的种群然后进入迭代循环。每次迭代对当前个体选中的 7 个站重新做 Voronoi 分区计算每个站覆盖的负荷再调用 VorCostCDEV 得到目标值。步骤调用函数输出读取候选点与负荷数据数据文件data 结构体生成初始种群main.mpop对当前解执行分区VoronoiT.m多边形顶点校核面积与负荷VoronoiArea.m各站服务负荷计算成本与网损VorCostCDEV.m个体目标值更新种群优化算法新一代个体% main.m 主循环示意 for iter 1:maxIter pop updatePopulation(pop); for i 1:popSize site data.site(find(pop(i,:) 1), :); [V, C] VoronoiT(site); for k 1:size(site, 1) verts V(C{k}, :); if ~any(isinf(verts(:))) areas(k) VoronoiArea(verts); end end cost(i) 0; for k 1:size(site, 1) loadsCell assignLoadsToCell(data.loads, verts); cost(i) cost(i) VorCostCDEV(site(k,1), site(k,2), ... Pcap(k), loadsCell, price, life); end end [bestCost, idx] min(cost); bestSite site(idx); end这段代码省略了 updatePopulation 的具体实现它可以是遗传算子也可以是粒子群速度更新。注意 VoronoiT 的输入只有当前选中的站点坐标而不是全部 34 个点否则未选中点也会划出多边形把 7 个站的服务区切碎结果毫无意义。assignLoadsToCell 的作用是判断负荷点属于哪个站可以用 matlab 的 inpolygon 实现也可以用 Voronoi 图自带的最近邻归属来判断。4. 从 34 个候选点选 7 个求解器、参数与收敛判断4.1 用 matlab 优化工具箱还是自己写粒子群34 选 7 的组合空间是 C(34,7)537 万看起来能穷举但每次评估都要跑一遍 Voronoi 分区和成本计算穷举并不现实。最省事的路径是直接用 matlab 优化工具箱自带的 ga 函数它支持整数约束可以处理 0-1 变量模型验证阶段足够。如果后续要把算法换成粒子群或模拟退火再自己实现也不迟。我在实际对比中发现ga 在 200 代以内基本能收敛到稳定解但偶尔会陷在局部最优。自己写粒子群时关键点是位置更新后需要把连续值映射回 0-1常见做法是加 sigmoid 函数并用 0.5 阈值判断建站与否同时用惩罚项保证正好选 7 个站。无论选哪种求解器VoronoiT、VoronoiArea、VorCostCDEV 这三个文件都不用改因为它们只负责评估一个给定方案的好坏。4.2 参数设置与约束处理用 ga 求解时代码框架可以写成nvars 34; lb zeros(1, nvars); ub ones(1, nvars); IntCon 1:nvars; % 所有变量都是整数配合 0/1 边界即为二进制 options optimoptions(ga, ... PopulationSize, 100, ... MaxGenerations, 200, ... CrossoverFraction, 0.8, ... EliteCount, 5, ... Display, iter); [x, fval] ga((x)CostWrapper(x, data), nvars, ... [], [], [], [], lb, ub, (x)conFun(x), IntCon, options);参数说明PopulationSize 是种群规模太小容易早熟太大会让每次迭代评估成本函数的次数暴涨MaxGenerations 是最大进化代数CrossoverFraction 控制交叉比例EliteCount 保证每代最优秀的个体不被破坏。conFun 里返回 ceq sum(x) - 7强制选 7 个站。0-1 变量不需要额外写约束边界 lb 和 ub 已经限死。参数建议值说明PopulationSize100组合规模不大100 足够MaxGenerations200观察 Display 输出决定是否增大CrossoverFraction0.8常用区间 0.7~0.9EliteCount5防止最优解被交叉变异破坏StallGenLimit5050 代无改善就停止CostWrapper 函数内部要做两件事一是调用 VoronoiT 等文件算出原始总成本二是对不满足约束的解施加惩罚。常见做法是 penalty 1e6 * abs(sum(x) - 7)直接把不符合 34 选 7 的解的目标值抬高这样优化器会自然避开。这里需要提醒惩罚系数不能设得过大否则可行域内外的目标差异被放大收敛曲线会非常难看也不能过小否则最终结果可能选出 6 个或 8 个站。另外ga 的初始种群是随机生成的如果不加 InitialPopulationMatrix 选项每次起点都不同。可以先用一个贪心算法生成初始解比如按负荷密度排序选前 7 个站然后塞进初始种群这样收敛速度会明显加快。4.3 排错结果不合理时先查这三处我拿到这类 matlab 程序时第一件事不是先看算法而是先跑一次默认参数确认 Voronoi 图能正常画出来。以下三个问题是初学者最容易碰到的。第一个问题是所有站聚在一堆或者部分站落在规划区域外。原因多半是缺少站间最小距离约束和边界约束。解决办法是在 CostWrapper 里加一个判断如果任意两个站距离小于 Dmin就返回一个很大的惩罚值让算法放弃这类方案。第二个问题是 Voronoi 面积出现 NaN导致目标函数全部变成 NaN。原因就是 3.1 节提到的 Inf 顶点没有裁剪。可以先用 matlab 画图验证voronoi(site(:,1), site(:,2))如果图上有延伸到无限远的射线说明需要做边界裁剪。裁剪可以用 inpolygon 逐点判断也可以把规划区域的四个角点加入 voronoin 的输入。第三个问题是下载的程序文件找不到函数。压缩包解压后没有添加到 matlab 路径或者当前文件夹不在工作目录都会报 Undefined function。解决方法是在 main.m 开头加一行addpath(genpath(pwd))或者手动右键文件夹添加到路径。这跟 matlab 版本无关纯路径问题但很多人会误以为是程序写错了。提示ga 或粒子群都是随机算法每次运行结果有波动是正常的。对比参数前先固定随机种子 rng(42)确保同一代码跑出同一结果否则你看到的差异是随机噪声而不是参数差异。5. 进阶把固定 7 个站改成动态定容并验证结果合理性5.1 把站点数量从固定值改成变量原模型固定 7 个站但实际规划需要回答“到底建几个站最划算”。最简单的改法是把约束从 sum(x)7 改成 sum(x)Kmax同时在目标函数里增加一项建站固定成本比如每个站加 100 万元固定投资。运行后如果某个站的收益盖不住固定成本优化器会自动放弃它。需要注意这时候 conFun 里的 ceq 变成了不等式约束ga 中需要写成 c sum(x) - Kmax并令 ceq []。5.2 画图验证 Voronoi 分区与负荷匹配选址结果出来以后不要只看一维成本数字建议用 matlab 画图检查用voronoi(site(:,1), site(:,2))叠加画出负荷点用不同颜色标记不同分区再在每个站旁标注容量。重点检查有没有多边形面积很大但容量很小的站这种站明显覆盖了超出能力的负荷说明容量约束写得不对或者 Voronoi 分区里混入了未选中的站点。验证代码可以这样写figure; voronoi(site(:,1), site(:,2)); hold on; scatter(data.loads(:,1), data.loads(:,2), 20, data.loads(:,3), filled); colorbar; for i 1:size(site,1) text(site(i,1), site(i,2), sprintf(station %d: %.0f kW, i, Pcap(i))); end这里用 scatter 的第三维给负荷点着色能直观看到充电需求密度高的区域是否被覆盖。如果某个颜色深的负荷簇恰好落在两个站的边界上就需要检查是不是站间距离约束太松调整 Dmin 后重新跑往往能满足该簇需求。5.3 扫描不同的 K 找成本拐点另一种验证方法是把 K 从 5 扫到 10分别记录最优目标值然后画总成本随 K 变化的曲线。如果 K7 时总成本明显低于 K6而 K8 比 K7 只低一点点说明 7 个站已经接近最优。把成本差除以新增站的年化投资就能判断多建一个站是否划算。扫描 K 时要注意每次运行都要固定 rng 种子否则不同 K 之间的随机噪声会盖过真实的成本变化。我一般会每个 K 跑 5 次取最小值再把最小值连成曲线这样得到的拐点更稳。碰到 U-shaped 曲线也不奇怪因为站太少网损高站太多投资高最低点就是规划期内的最佳站数。本文还有配套的精品资源点击获取