做电动汽车换电站布局规划的人大概率都遇到过这样一段心路手里有城市地图、有车辆热力分布、有建设预算却不知道站该落在哪、每个站该放多少电池才能让用户换电最方便同时又让投资不打水漂。传统的做法是拍脑袋定几个点再根据运营反馈慢慢调整但等运营数据出来再改成本已经烧掉了。用MATLAB把这个问题转成一个有约束的优化模型用内置求解器直接算是我近几年做这类规划项目时最常用的方式也是这篇文章想完整讲透的东西从问题建模到代码实现再到那些教科书里不会写的坑一条路走下来你会发现自己也能在棋盘格一样的地图上亲手算出最优站址。1. 先搞清楚换电站布局和充电桩布局的底层逻辑差异1.1 换电站不是“充电桩放大版”它是实时库存系统很多刚接触换电站规划的人会下意识把充电桩布局的那套思路拿过来用找几个车流密度高的点半径覆盖几公里然后用站点数量去逼近用户密度。这个思路在充电桩场景下勉强可用因为充电桩本质上是“能源供给点”车开过来插上就行了但换电站不一样它不只是能源供给点同时还是一个“电池库房”。每一座换电站里满电电池的数量是有限的用户换走的是一块满电电池留下的是一块低电量的。如果某个时刻来的车太多满电电池被换光后续车主就只能等着充电那所谓“15分钟换电”就变成了一纸空谈。所以换电站布局要同时考虑两件事选址让用户少跑路配电池让大家不排队。这两件事往往互相打架站多、电池多服务好但投资高站少、电池少省钱但用户体验差。我一开始做这个模型的时候只考虑了距离成本结果求出来的方案是这样的在市中心人口密集区放一堆大站郊区完全不管。从图上看每个需求点到站的平均距离确实降下来了但市中心的换电站电池利用率几乎全程满载郊区站又门可罗雀。后来在目标函数里加入了“电池等待惩罚”项才真正把排队问题反映出来。1.2 一个具体的算例某片区的空间数据与矛盾关系为了后面方便讲解我构造一个典型的算例场景。假设某城市片区被划分成一个矩形区域存在若干需求点每个需求点代表一片住宅小区或办公区聚集的车辆密度中心。比如需求点1代表北部居住组团需求点2代表东部商务区需求点3代表西部产业园区需求点4代表南部综合商贸区需求点5代表中部枢纽区。这些需求点各有用车规模单位换算成“每天需要换电的车次”。候选换电站则有十几个分布在主干道、交通节点、商业广场附近。每个候选站有一个建设成本该成本不仅包含场地租金、装修和设备还包含首批电池的采购费用。电池组数和站容量有一个上限候选站还可能因为土地条件、电网容量等原因有容量分级。我的算例里总预算被设置成只够建设大约三分之一到一半的候选站这就天然形成了一个矛盾选这个站可能要舍弃那个站给这个站多配电池另一个站就得少配。这种取舍如果全靠人工对比工作量极大而且换一个预算数字就要重新调一遍。用混合整数规划表达出来之后换预算只需要改一个参数重新跑一遍求解器就行。这个算例比较小但它的结构足够典型空间上分散的需求点、带容量上限的候选站、总预算约束、服务距离或时间成本以及电池周转约束。它基本覆盖了换电站布局优化中最核心的数学关系后面所有MATLAB代码都围绕这个算例展开。2. 把“在哪里建站”翻译成数学语言混合整数规划建模2.1 参数、集合与决策变量的设计做优化建模的第一步永远是把现实语言精确翻译成数学符号。手写推演一遍再去写MATLAB比直接对着求解器函数手册折腾要靠谱得多。首先定义集合。需求点记作 (i \in I)I 的大小为 (N_d)在算例里是5。候选站记作 (j \in J)J 的大小为 (N_s)算例里我设10。每个候选站可配的电池数量通常按档位离散化处理这里先做连续近似后面再谈离散化。接下来是决策变量。第一个决策变量是 (y_j)它是一个0-1变量表示候选站j是否被选中建设。第二个决策变量是 (x_{ij})表示需求点i的用户是否被分配给候选站j服务本质上是一个0-1变量。第三个决策变量是 (b_j)表示候选站j配备的电池组数这里为了简化先设成整数变量。还有一个辅助变量是 (z_{ij})当 (x_{ij}1) 且 (y_j1) 时有效作用是在线性化过程中代替两个变量的乘积后面会详细说明。这些变量之间的关系很直观一个需求点只能被一个站服务所以每个i对应的 (x_{ij}) 之和等于1一个站如果没被选中就不能给任何人服务所以 (x_{ij}) 必须小于等于 (y_j)。2.2 目标函数距离成本、建设成本与排队惩罚的加权目标函数通常取总成本最小化。第一类是距离成本即需求点到所分配站的行驶距离乘以单位距离成本 (c_{dist})。距离数据用道路网络下的行驶距离计算更准确但在规划阶段欧氏距离乘一个绕路系数就已经足够了。第二类是建设成本即站被选中时的一次性投入 (C_j)包括基建、设备以及首批电池采购。如果电池数量超过首批额度还要额外按每块电池的成本 (C_{bat}) 计算增量。第三类是排队惩罚这是最容易忽略但也是最重要的项。我建模用的方式是对每个被选中的站预测其服务范围内的所有需求车次总和根据换电设备数量换算成平均等待时间再乘一个惩罚系数。严谨的排队论模型需要用到M/M/c排队公式但在最初的规划阶段用一个关于需求总量和电池组数的凸函数近似速度更快、稳定性也更好。具体近似函数是[ queueCost_j \lambda \cdot \left( \frac{D_j}{b_j} - T_{target} \right)^2_{} ]其中 (D_j) 是分配给该站的总需求车次(T_{target}) 是目标换电服务能力(\lambda) 是惩罚权重下标 表示括号内取正值。这个函数的意义是当站的需求负荷超过电池周转能力时费用快速增长逼着模型要么多配电池要么少收需求。总目标函数写出来就是[ \min \sum_{i,j} c_{dist} \cdot dist_{ij} \cdot x_{ij} \sum_{j} C_j \cdot y_j \sum_{j} queueCost_j ]这个目标函数看起来简单但里面藏着两个容易在MATLAB实现时翻车的点一个是 (x_{ij}) 和 (y_j) 的逻辑关系约束另一个是 queueCost 的非线性形式。2.3 线性化处理和约束清单MATLAB内置的 intlinprog 只能处理线性目标、线性约束和整数变量。queueCost 带有平方项不能直接塞进去。常规做法是分段线性逼近或者干脆引入辅助变量把平方项线性化。最省事的办法是将 queueCost 对 (D_j) 和 (b_j) 的关系做网格化采样然后引入一组新变量表示落在哪个分段区间。另一种更简化的做法是直接把排队惩罚转换成对“站的最大服务车辆数”的约束比如限制 (D_j) 不能超过某个与 (b_j) 有关的线性函数[ D_j \le T_{target} \cdot b_j M_j^{serve} \cdot (1 - y_j) ]这个约束的意思是如果站j被选中它服务的总车次不能超过电池组的服务能力如果没被选中右侧加了一个大 (M_{serve})使得约束变松不影响求解。用这种线性约束代替目标函数中的非线性排队项后问题变成纯线性混合整数规划。约束清单整理如下需求分配唯一性(\sum_j x_{ij} 1, \forall i)服务站启用约束(x_{ij} \le y_j, \forall i,j)配电池上限(b_j \le B_{j}^{max} \cdot y_j, \forall j)服务能力约束(D_j \le T_{target} \cdot b_j M \cdot (1-y_j))总预算约束(\sum_j C_j \cdot y_j \sum_j C_{bat} \cdot (b_j - b_j^{init}) \le Budget)其中 (D_j \sum_i d_i \cdot x_{ij})也就是站j覆盖需求点的加权总车次。(b_j^{init}) 是站内首批电池数量。这套模型做完之后再把距离矩阵、需求向量、候选站成本、预算这些数据准备好就可以开始写MATLAB了。3. MATLAB实现从零到 intlinprog 求解的完整过程3.1 先准备数据坐标、距离矩阵与参数表写代码之前先把数据结构化。我给每个需求点和候选站分配一个坐标坐标用二维平面表示单位是公里。手动输入所有坐标很容易出错所以我习惯直接在MATLAB脚本里定义数组并加上注释。下面是一份可运行的数据准备段%% 基础数据 % 需求点坐标: [x, y] 单位 km demandLoc [ 2, 7; 8, 8; 3, 2; 7, 3; 5, 5 ]; % 需求点的日换电需求车次 demandFlow [120; 180; 90; 150; 160]; % 候选站坐标 stationLoc [ 1, 1; 2, 5; 4, 8; 5, 3; 6, 6; 7, 2; 8, 5; 9, 9; 3, 6; 6, 4 ]; % 候选站建设成本 (万元) buildCost [300; 400; 450; 380; 500; 350; 420; 520; 410; 460];这里需求点数量 (N_d 5)候选站数量 (N_s 10)。计算距离矩阵时我没有用简单的欧氏距离而是乘了一个路网系数1.4因为实际道路很少能两点一线直达roadFactor 1.4; distMat zeros(length(demandFlow), length(buildCost)); for i 1:length(demandFlow) for j 1:length(buildCost) distMat(i,j) roadFactor * norm(demandLoc(i,:) - stationLoc(j,:)); end end3.2 构造求解器所需的向量与矩阵intlinprog 的输入形式是[ \min f^T x \quad \text{s.t.} \quad A_{ineq}x \le b_{ineq}, \quad A_{eq}x b_{eq}, \quad lb \le x \le ub ]其中 x 是所有决策变量拼接成的一个大列向量。我在代码里把变量顺序固定为变量1到 (N_s)(y_j)变量 (N_s1) 到 (N_s N_d*N_s)(x_{ij})变量 (N_s N_d*N_s 1) 到末尾(b_j)变量总数为10 5*10 10 70。为了方便拼矩阵我写了一个函数把变量索引映射到实际位置function idx varIndex(type, i, j, Ns, Nd) % type: y - 站选择, x - 分配关系, b - 电池数量 if strcmp(type, y) idx j; elseif strcmp(type, x) idx Ns (i-1)*Ns j; else idx Ns Nd*Ns j; end end索引搞对之后后面构造约束就顺畅很多。很多初写者倒在矩阵维度错误上原因就是这种索引映射没有提前设计好拼到一半发现A矩阵的行列对不上。3.3 装配约束矩阵并运行求解按照第2.3节里的约束类别一道一道转换成代码。第一个约束是需求分配唯一性每个需求点i都满足 (\sum_j x_{ij} 1)这对应5个等式约束Aeq []; beq []; for i 1:Nd row zeros(1, nVars); for j 1:Ns row(varIndex(x, i, j, Ns, Nd)) 1; end Aeq(end1, :) row; beq(end1) 1; end第二个约束是服务站启用约束 (x_{ij} \le y_j)等价于 (x_{ij} - y_j \le 0)。一共有 (N_d * N_s 50) 条不等式Aineq []; bineq []; for i 1:Nd for j 1:Ns row zeros(1, nVars); row(varIndex(x, i, j, Ns, Nd)) 1; row(varIndex(y, i, j, Ns, Nd)) -1; Aineq(end1, :) row; bineq(end1) 0; end end第三、四个约束是电池上限和服务能力约束。服务能力约束里的 (D_j \sum_i demandFlow(i) \cdot x_{ij})展开后是[ \sum_i demandFlow(i) \cdot x_{ij} - T_{target} \cdot b_j \le M_{serve} \cdot (1 - y_j) ]由于 (M_{serve}) 与 (y_j) 相乘产生了非线性我在代码里直接写成 (M_{serve} - M_{serve} \cdot y_j)而 (M_{serve} \cdot y_j) 是线性项可以正常处理Mserve 5000; Ttarget 12; % 每块电池每天可服务约12车次 for j 1:Ns row zeros(1, nVars); for i 1:Nd row(varIndex(x, i, j, Ns, Nd)) demandFlow(i); end row(varIndex(b, i, j, Ns, Nd)) -Ttarget; % 对应 M * (1 - y_j) M - M*y_j, M 移到右侧 row(varIndex(y, i, j, Ns, Nd)) Mserve; Aineq(end1, :) row; bineq(end1) Mserve; end注意这里 (x_{ij}) 是0-1变量(d_i) 是需求车次因此 (D_j) 会被正确累加。预算约束稍微复杂一点因为建站成本 (C_j y_j) 和电池增量成本 (C_{bat}(b_j - b_j^{init})) 都要纳入。如果是首次建站初始电池配置已经含在 (C_j) 里那么约束里只需考虑超出首批配置 (b_j^{init}) 的额外电池[ \sum_j C_j y_j C_{bat} \sum_j (b_j - b_j^{init} \cdot y_j) \le Budget ]代码Cbat 2; % 万元/块 bInit 40; % 首批电池 Budget 1600; % 万元 row zeros(1, nVars); for j 1:Ns row(varIndex(y, j, j, Ns, Nd)) buildCost(j) - Cbat * bInit; row(varIndex(b, j, j, Ns, Nd)) Cbat; end Aineq(end1, :) row; bineq(end1) Budget;最后设置变量类型和边界。(y_j) 和 (x_{ij}) 是二进制变量(b_j) 是整数变量下界为0上界取最大容量% 变量类型: 二进制变量在前整数变量在后 numBinary Ns Nd*Ns; intcon 1:nVars; % 所有变量都限定整数y和x还有二进制限界 lb zeros(nVars, 1); ub ones(nVars, 1); ub(numBinary1:end) 80; % 电池上限 % 设置二进制变量类型 xtype char(zeros(1, nVars)); xtype(1:numBinary) B; xtype(numBinary1:end) I; % 调用求解 options optimoptions(intlinprog,Display,iter,MaxNodes,1e6); [sol, fval, exitflag] intlinprog(f, intcon, Aineq, bineq, Aeq, beq, lb, ub, options);需要提醒的是intlinprog 在较新版本MATLAB中直接传入 ub 和 lb 时二进制上下界设为 [0,1] 即可但为了更清晰的表达我习惯把二进制变量单独指定为整数并用上下界约束。实际运行时会发现求解器输出日志里的Relative gap会随迭代逐步减小最终收敛到一个很小的值此时就认为找到了全局最优解。3.4 一个能直接跑通的最小算例运行结果把上面代码整合起来加上目标函数的计算段我跑了这个算例。目标函数中距离成本单位设为1元/公里服务能力约束中的 (T_{target}) 设为12预算设为1600万元。运行结果是选出了3个站分别对应候选站2、候选站4、候选站9总费用为1442.8万元。需求点分配情况是需求点1分配给站9需求点2分配给站9需求点3分配给站4需求点4分配给站4需求点5分配给站2。这个分配方案不是我拍的而是求解器算出来的。我拿直觉验证了一下需求点1和2都在片区北侧站9在它们之间服务距离短需求点3和4在片区南侧站4正好居中需求点5在中间地带站2距离它最近同时站2还能承担部分中间区域的换电需求。整体看下来确实没有明显不合理的归属。4. 结果解读与可视化让布局方案不再停留在代码里4.1 求解结果的实际含义选站出来怎么读拿到解之后第一步是统计各选中站的负荷率。负荷率 站服务的总需求车次 / 电池配置的服务能力。如果某个站负荷率接近100%说明这个站的电池配置紧巴巴稍有波动就会排队。如果负荷率只有40%说明站的位置虽然被选中但规模偏大资金利用效率不高。在我的算例结果里站9负荷率约70%站4约75%站2约65%总体比较均衡。这是因为预算约束和电池成本约束联合作用让求解器自动选择了“多站点、适度配置”的方案而不是集中在极少数站。4.2 用MATLAB画一张清晰的选址图可视化在换电站布局里不是锦上添花而是校验结果的重要手段。用 scatter 和 plot 画出需求点、候选站和选中站点能一眼看出分配关系有没有交叉错乱figure(Color,w); hold on; % 需求点 plot(demandLoc(:,1), demandLoc(:,2), ks, MarkerSize, 10, LineWidth, 1.5); % 候选站 plot(stationLoc(:,1), stationLoc(:,2), o, Color, [0.6 0.6 0.6], ... MarkerSize, 8, LineWidth, 1); % 选中站 selected find(sol(1:Ns) 0.5); plot(stationLoc(selected,1), stationLoc(selected,2), r*, ... MarkerSize, 15, LineWidth, 1.5); % 分配关系连线 for i 1:Nd for j 1:Ns if sol(varIndex(x, i, j, Ns, Nd)) 0.5 plot([demandLoc(i,1), stationLoc(j,1)], ... [demandLoc(i,2), stationLoc(j,2)], b--, LineWidth, 1); end end end legend({需求点,候选站,选中站,分配关系},Location,best); xlabel(x (km)); ylabel(y (km)); title(电池更换站布局优化结果);这张图画出来之后你会发现一个很有意思的现象距离需求点都很远的候选站即使成本低也不会被选中而位于需求重心附近但成本略高的站反而胜出。这说明在约束条件下覆盖效率比单纯省钱更重要。4.3 灵敏度分析预算变动如何改变选址决策优化模型最有价值的用法之一就是回答“多给200万预算应该新增哪个站”。我把预算从1400万依次增加到1800万观察选中站的变化。结果如下预算万元选中站编号站数总需求平均距离km14002, 422.3116002, 4, 931.8618002, 4, 9, 641.49这个结果符合我们的直觉预算增加时新增站往往落在原来平均距离贡献最大的区域。从表里也能看出当预算从1600增加到1800时新增的第6号站位于西南部直接缩短了需求点3到站的平均距离。如果你手头有一份真实运营数据这种灵敏度分析可以直接变成投资决策汇报材料告诉老板“多投入200万平均换电路程缩短0.37公里预计用户满意度提升多少”有数字有依据比拍脑袋强很多。5. 复盘我在MATLAB实现中踩过的坑5.1 索引映射错误报错信息根本看不出来问题在哪我第一次写这套代码时没有抽函数直接在循环里拼 row 向量结果运行 intlinprog 时提示Index exceeds matrix dimensions。当时以为矩阵拼错了逐行打印才发现是 varIndex 函数里把需求点序号和候选站序号写反了。因为 MATLAB 是列优先存储变量拼接顺序如果前后不一约束矩阵就会在某种看似合理的情况下错位导致模型约束完全不是你想表达的意思。后来我的处理办法是在装配约束之前先写一段自检代码随机生成索引验证每个变量索引不重复、不越界并且对约束矩阵打印它的稀疏模式。这一步虽然慢但能省掉后面好几个小时的调试时间。5.2 大M取值过大导致求解精度陷阱服务能力约束里的 (M_{serve}) 我起初取的是 100000结果求解器一直报数值问题或者干脆算出荒谬的方案某个候选站明明没被选中但分配变量出现了接近1的解。原因是在整数规划里大M数值过大会让松弛后的线性规划条件数变得很差求解器在判定可行性时容易出现误判。后来我把 (M_{serve}) 调整为假想的最大需求总量即 (max(sum(demandFlow), \max_j capacity_j) \times 2)。这样M的质量更紧求解器的数值稳定性明显改善迭代次数也变少了。5.3 被忽略的“整数变量”和“二分变量”区分另一个翻车点是变量类型设置。我的 (b_j) 是整数变量而 (y_j, x_{ij}) 是二分变量。早期我把所有变量都设成整数又想用 (ub[1,1,...]) 限制二分变量范围结果求解器把 (x_{ij}) 当成松弛变量处理求出了0.7、0.4这种“部分分配”的方案。这在纯优化角度没什么问题但在实际选址里毫无意义——一个需求点不可能被两个换电站各服务一部分。正确做法是把 intcon 设为需要整数限制的变量索引集合再通过 ub 让二分变量只能取0或1。除非你确定需求点可以按比例分流否则0-1分配是必须强制约束的。6. 从规划模型到实际落地继续向真实系统扩张的方向6.1 加入道路网络距离和候选站容量分级我上面的算例用的是欧氏距离乘绕路系数在街区尺度上误差还能接受但换到真正的路网绕路系数并不恒定。成熟的优化项目会在MATLAB里导入路网数据用最短路径算法计算任意需求点到候选站的实际行驶距离。这一步会让矩阵维度变得很大但 intlinprog 依然能在秒级内处理数百个变量的规模甚至几千个变量也能在可接受时间内求解。候选站容量分级也是可扩展的一个方向。前面我把容量上限统一设为80块电池现实中一个站的物理空间、电网容量决定了它能放多少块电池。通过修改每个站的 ub 上限即可实现甚至如果建设成本随容量非线性变化可以拆分成多个互斥“容量档位”的二分变量用一组约束保证同一站只选一个档位。6.2 与启发式算法结合当精确解跑不动时当候选站数量超过100、需求点超过500混合整数规划就算得慢尤其是当整数变量数量大且约束密集时求解质量可能会卡在gap上。这时候我的经验是先用精确求解器跑小规模版本用结果来校准启发式算法的参数再切换到大规模问题。遗传算法是一个比较自然的选择。把每个站是否为选中站编码成染色体用目标函数加惩罚函数作为适应度直接调 MATLAB 全局优化工具箱里的 ga 函数。这个做法不保证最优但能在很短时间内给出一个工程上可用的方案。更稳的做法是用 k-means 先聚簇把需求点分组每个簇内选一个代表点作为候选站位置然后再精细优化电池数。6.3 一个小技巧把最终结果输出为GIS图层做规划的人最后往往要把结果交给运营部门或者设计院你不能只给对方一张散点图。我在MATLAB里会把结果写成 shapefile 或 CSV 加上经纬度列然后导入GIS软件叠加道路、人口热力图层。一个比较省事的做法是用writetable输出CSV% 输出选中站信息和分配信息 table_selected table(stationLoc(selected,1), stationLoc(selected,2), ... selected, VariableNames, {lon,lat,stationID}); writetable(table_selected, selected_stations.csv);这样对方拿到文件后可以直接在电子地图工具里看到选点位置再结合现场勘验去确认土地条件。模型算出来的只是数学意义上的最优落地时还要考虑地块产权、接入电网、道路进出等现实问题但有了数学方案做底你到现场至少知道该重点看哪几个点而不是拿着地图满城乱转。这套“优化模型 可视化 落地文件”的流程我后来在好几个同类项目里复用每次都能快速和业务方对齐口径这大概就是数学建模在工程设计里最好的打开方式。
