简介：面向电动汽车充电站规划场景的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这段代码把决策变量 x（34 维 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 | 候选点数量 | 本例 N=34 |
| K | 计划建站数 | 本例 K=7 |
| x | 34 维 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_loss，main.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.m | pop |
| 对当前解执行分区 | 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 已经限死。
| 参数 | 建议值 | 说明 |
|---|---|---|
| PopulationSize | 100 | 组合规模不大，100 足够 |
| MaxGenerations | 200 | 观察 Display 输出决定是否增大 |
| CrossoverFraction | 0.8 | 常用区间 0.7~0.9 |
| EliteCount | 5 | 防止最优解被交叉变异破坏 |
| StallGenLimit | 50 | 50 代无改善就停止 |
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 变化的曲线。如果 K=7 时总成本明显低于 K=6，而 K=8 比 K=7 只低一点点，说明 7 个站已经接近最优。把成本差除以新增站的年化投资，就能判断多建一个站是否划算。
扫描 K 时要注意，每次运行都要固定 rng 种子，否则不同 K 之间的随机噪声会盖过真实的成本变化。我一般会每个 K 跑 5 次取最小值，再把最小值连成曲线，这样得到的拐点更稳。碰到 U-shaped 曲线也不奇怪，因为站太少网损高，站太多投资高，最低点就是规划期内的最佳站数。
本文还有配套的精品资源，点击获取
