考虑不确定性的含集群电动汽车微电网随机优化调度实现-编程实验室

做微电网调度的朋友应该都有同感：这两年随着电动汽车保有量持续上升，原本“电源—储能—负荷”三段式的微电网优化调度问题，硬生生多了一个变量维度。车不再是单纯的负荷，它还自带电池，随时可能作为移动储能参与充放电，这就是课题里“集群电动汽车”的由来。再加上风电、光伏、用户负荷本身就飘忽不定，调度决策如果不把这些不确定性考虑进去，得到的结果往往只能停留在论文图表里，到了现场执行很容易出现功率缺额、弃风弃光甚至频率波动的问题。

这篇博文要聊的，就是“考虑不确定性的含集群电动汽车并网型微电网随机优化调度”这个课题的完整实现路径。包括集群电动汽车怎么建模、风光伏和负荷的不确定性如何用场景刻画、随机优化模型怎么搭建，以及Matlab代码如何一步步落地求解。我会尽量用做工程的口吻讲清楚每一个选择背后的原因，也会把调试中踩过的坑一并写出来。不管你是正在做微电网相关毕设、准备入门随机优化，还是已经写了代码但结果总是不对，这篇文章都值得花时间看完。

1. 先弄清楚：这个课题到底在解决什么问题

1.1 集群电动汽车为什么让调度变得棘手

先铺垫一下传统微电网调度是什么。并网型微电网里的主要电源是风电、光伏、微型燃气轮机，储能负责平抑波动，负荷端是普通电负荷和可调负荷。调度问题本质上是一个多时段、多单元的能量管理问题：在满足各种运行约束的前提下，以最小成本或最大收益决定每个单元每个时段的出力。这个框架在“没有EV”的时代已经非常成熟，主流做法是确定性优化加滚动修正，工程上跑得很顺。

但电动汽车一进来，局面就变了。EV对调度系统来说非常特别：它既是负荷（充电时），又是电源（V2G放电时），还是储能设备（停在原地时）。更棘手的是它的“身份”会随车主行为实时变化——什么时候接入电网、什么时候拔枪离开、初始电量是多少、今天要跑多少里程，这些都带随机性。很多初学者在第一版建模时习惯把EV当作普通储能直接写进模型，结果调度出来的方案要求某个时段有大量EV放电，但实际上那个时段大部分车都在路上跑，方案根本无法执行。这个问题做仿真时不容易暴露，一到实际项目就原形毕露。

拿我自己做过的案例来说，某居民区微电网假设有200辆可调度EV，晚高峰18:00到21:00之间真实“在网且可调度”的EV可能只有100辆出头，而且它们的目标SOC各不相同。如果不把“接入时间、离开时间、出行需求”这三件事建模清楚，聚合出来的可调度容量可能比真实值高出40%到60%，调度结果的安全裕度会被严重高估。这也是为什么所有正经的EV集群研究，第一步都必须做“个体行为建模+聚合”。

1.2 四类不确定性：不处理就等着“翻车”

不确定性是这个课题的核心关键词，它主要来自四个地方：风电出力的随机性（风速随机且预测误差随时效增大）、光伏出力的随机性（云层遮挡导致出力波动剧烈）、负荷的随机性（用户用电行为波动）、以及EV行为随机性。如果微电网参与电力市场交易，还要算上实时电价的不确定性。我把常用建模方式整理成一个表，方便大家对照：

不确定性来源	典型概率分布	说明
风电出力	Weibull分布描述风速，再经功率曲线转换	短期预测误差常近似为正态分布
光伏出力	Beta分布描述光照强度	同样常用预测误差正态近似
负荷	正态分布	预测误差随预测周期增大而增大
EV群体行为	到达/离开时间用分段正态或泊松分布，日行驶里程用对数正态分布	初始SOC可由行驶里程反推
实时电价	正态分布或场景集	参与市场时可通过场景法处理

为什么不能对不确定性视而不见？做个简单估算：假设并网微电网的联络线功率上限是500kW，某时段风光预测出力600kW，若按确定性思想调度，这个时段可以把燃气轮机停机，让负荷全部由风光供。但实际风光出力可能只有450kW，这时为了维持功率平衡，必须立刻从电网多买电，购电成本可能瞬间拉高；如果联络线已经满载，还会被迫切负荷。一个成熟的随机优化模型，就是通过场景把这种“差一点就翻车”的情况提前纳入优化，在成本和风险之间找平衡点。

2. 核心建模：把实际问题变成可求解的数学问题

2.1 集群电动汽车：聚合模型是关键

单体EV建模并不复杂，接入电网后可以看作一台小型储能，关键参数有：电池容量（一般20到100kWh）、最大充放电功率（交流慢充3.3到7kW，直流快充50到120kW）、充放电效率（0.9到0.95）、初始SOC、离开时目标SOC。车主的出行需求通过“目标SOC”体现：当天行驶里程越长，接入时需要补充的电量就越多。

集群聚合的目标，是把几百辆单车用一个等值参数描述，但必须保留“时间段”的信息。假设调度周期为24个时段、每时段1小时，定义t时段在网的EV集合为Ω(t)，那么聚合充电功率上下限就不是一个常数，而是由每个时段在网车辆数决定的时变区间。以充电方向G2V为例，聚合模型可以写成：

聚合功率下限：P_ev_min(t) = Σ P_charge_i_min(t)，i ∈ Ω(t)
聚合功率上限：P_ev_max(t) = Σ P_charge_i_max(t)，i ∈ Ω(t)
能量递推：E_ev(t) = E_ev(t-1) + (η_c * P_ev_charge(t) - P_ev_discharge(t) / η_d) * Δt
能量上下限：E_ev_min(t) ≤ E_ev(t) ≤ E_ev_max(t)
调度期末约束：E_ev(T_end) ≥ Σ E_target_i，保证每辆车离开时达到目标SOC

如果你用确定性调度，直接把聚合参数当常数填进去；如果做随机规划，就要先对EV群体里的随机个体行为做蒙特卡洛模拟，再聚合出每个场景下的P_ev_min、P_ev_max、E_ev_min、E_ev_max等参数。我建议初学者先把单车的“时间窗约束”想清楚，再动手写聚合，否则很容易写出“每辆车全天都能充放电”的错误假设。

另外，EV集群并不一定全部可调度。实际中，很大一部分车主是“即插即充”模式，只有少部分愿意参与V2G，所以模型里通常还要加一个“可调度比例”参数，比如只有30%的车辆参与优化调度，其余按固定充电曲线处理。这个参数对充电桩运营商的实际运营至关重要，也是审稿人或答辩老师很喜欢追问的细节。

2.2 风、光、负荷的不确定性建模

不确定性建模理论上分两层：一层是物理过程产生的不确定性，另一层是预测误差导致的不确定性。调度模型里通常用的是“预测误差”那一层，因为调度面向未来，我们手里只有预测曲线和误差分布。

风电方面，风速可以用两参数Weibull分布描述，Matlab里wblpdf、wbllcdf可以直接调用，不需要手写公式。风速通过风机功率曲线映射为出力，切入风速一般在3到4m/s，额定风速11到13m/s，切出风速22到25m/s。场景生成时直接对风速做蒙特卡洛采样，再经功率曲线转到出力即可。

光伏出力常用Beta分布描述太阳辐照度的随机性，再由辐照度和温度通过光伏模型得到出力。不过在调度类代码里，大家更习惯的做法是假设预测误差服从均值为0的正态分布，把“预测值+误差”当作出力场景。例如某时段光伏预测出力是P_pv_forecast，误差标准差取预测值的10%到20%，采样时生成P_pv_scenario = P_pv_forecast + error。负荷的不确定性与光伏类似，常用正态误差模型，标准差按负荷规模取2%到5%。这个数值来自工程经验，没有统一标准，需要根据实际历史数据确定。如果你手头有历史出力数据，建议先用histfit看一眼分布形态，再决定用正态分布还是厚尾分布。

到这里模型基础已经搭完：EV聚合给出“可调度资源”的时变参数，风光负荷给出随机出力的分布描述。下一步就是把随机性放进优化框架。

3. 随机优化框架：用场景集逼近真实世界

3.1 为什么选随机优化，而不是鲁棒优化或确定性优化

调度领域处理不确定性有三种主流思路：确定性优化、鲁棒优化、随机优化。用表格对比一下更直观：

方法	不确定参数处理	目标	优点	缺点
确定性优化	用预测值代替	单情景成本最低	简单、求解快	容易低估风险
鲁棒优化	区间/盒式集合	最坏情况成本最低	结果保守可靠	往往过于保守，经济性差
随机优化	场景集	期望成本最低	经济性与可行性平衡	需要合理的场景概率模型
分布鲁棒优化	模糊集	最坏分布下的期望最优	兼顾鲁棒与随机	模型复杂、求解难度大

选择随机优化最核心的原因是工程语义明确：EV到达时间、风机出力这类不确定性可以通过历史数据拟合出概率分布，恰好适配场景法。随机优化的目标函数不是最小化某一个场景下的成本，而是最小化所有场景下运行成本的数学期望，同时在每个场景下都要保证功率平衡和运行约束成立。这比鲁棒优化更“会过日子”——它不为“几乎不可能发生的最坏情况”全额买单，但通过场景概率加权，给大概率事件足够权重。

举一个电动车的例子。鲁棒优化会为“晚上20:00全部200辆车同时充电”这种极端场景设计调度方案，为了避免它，整个白天可能都在压低储能放电、维持高热备，成本很高；而随机优化看到的是“20:00大概率有120辆车在网、其中30%愿意参与调度”，据此做的日常决策更经济。当然这也带来一个要求：场景集必须能代表真实分布，否则随机优化同样会给出“过于乐观”的方案。

3.2 场景生成与削减：计算量和精度之间的平衡

蒙特卡洛采样生成场景是随机优化的第一步。假设风速、光照、负荷、EV行为都有各自的概率模型，那么一次采样就生成一整天的时序曲线，构成一个场景。常规做法是先采样几百到几千个原始场景（比如1000个），然后做场景削减，把相似场景合并，最终保留20到50个代表性场景，并给每个场景赋予概率。

场景削减最常用的是同步回代消除法，思路概括为：

计算所有场景两两之间的欧氏距离或Kantorovich距离；
找到“与其他场景距离最近”的场景，如果它和最近的另一个场景距离最小，就把其中一个删掉；
被删场景的概率加到离它最近的那个场景上；
重复直到剩余场景数到达预设值。

Matlab里可以手写这个函数，也可以用现成代码包。很多论文用kmeans做聚类削减，效果也不错，但kmeans的缺点在于聚类中心不一定是真实场景点，可能产生“物理上不可能”的代表场景，比如某个时段出力超过装机容量。同步回代法没有这个毛病，所以更推荐。

实际采样代码就几行：

nsample = 500; wind_scenarios = P_wind_forecast + wind_std * randn(nT, nsample); pv_scenarios = P_pv_forecast + pv_std * randn(nT, nsample); load_scenarios = P_load_forecast + load_std * randn(nT, nsample);

削减多少场景合适？我做过一组测试：同样的EV聚合数据和预测曲线，原始场景500个，削减到10、20、30、50分别求解，结果发现20个以上时优化目标已经基本稳定，10个以下开始出现明显偏差。综合计算量考虑，日常工程取20到30个场景是比较平衡的区间。注意每个场景在目标函数里都有概率权重，所以削减后别忘了给每个保留场景赋归一化概率，否则期望成本计算就是错的。

4. Matlab工程实现：从公式到可运行代码

4.1 代码结构和数据准备

我不太推荐把所有代码写在一个脚本里，虽然看起来“一键运行”，但调试起来后悔莫及。建议按模块拆文件：

main.m：主程序，负责加载数据、调用场景生成和优化函数、汇总结果和画图；
data_generation.m：生成风电、光伏、负荷预测曲线和EV集群参数，输出24小时时序数据；
scenario_generation.m：根据概率模型生成原始场景；
scenario_reduction.m：实现同步回代消除法，返回削减后的场景及概率；
ev_aggregation.m：给定单车EV参数，计算各时段聚合可用功率上下限和能量上下限；
optimization_model.m：用YALMIP构建随机优化模型并调用求解器；
plot_results.m：绘制机组出力、EV充放电、储能SOC、联络线功率等结果图。

数据准备阶段有两个容易踩坑的点。第一，时间索引要统一，所有数组都应该是24×1或1×24，别让光伏、负荷、EV参数之间差一个索引，否则约束写错时极难排查。第二，基准单位全部统一，比如功率用kW，储能容量用kWh，电价用元/kWh，单位混用会引起数量级错误，而且这类错误在优化问题里不会直接报错，只会让结果安静地“离谱”。

4.2 目标函数与约束的YALMIP实现

在随机优化里，需要区分两大类决策变量：第一阶段决策是现在必须确定、不随场景变化的变量，通常是和电网达成的购售电计划、机组启停状态；第二阶段决策是随场景变化的变量，指每个场景下各设备的具体出力、EV充放电功率等。两阶段随机规划的典型形式是：

min c' * x + Σ_s p_s * Q_s(x, ξ_s)

其中Q_s(x, ξ_s)是场景s下的运行成本函数。在YALMIP里定义变量时就要区分场景下标，核心代码大致长这样：

% 场景数据 ns = 20; % 削减后的场景数 p = redu_prob; % 场景概率，1×ns % 第一阶段决策变量 x_unit = binvar(nT, nG); % 机组启停状态 % 第二阶段决策变量，带场景下标 P_g = sdpvar(nT, nG, ns, 'full'); % 燃气轮机出力 P_bat = sdpvar(nT, 1, ns, 'full'); % 储能功率 E_bat = sdpvar(nT, 1, ns, 'full'); % 储能SOC P_ev = sdpvar(nT, 1, ns, 'full'); % EV聚合功率，正为充电 P_buy = sdpvar(nT, 1, ns, 'full'); % 向电网购电 P_sell = sdpvar(nT, 1, ns, 'full'); % 向电网售电 P_curtail = sdpvar(nT, 1, ns, 'full'); % 弃风弃光量 % 目标函数：期望运行成本 objective = 0; for s = 1:ns objective = objective + p(s) * ( ... sum(sum(c_g * P_g(:,:,s))) + ... % 燃气机组燃料成本 sum(c_buy * P_buy(:,s)) - ... % 购电成本 sum(c_sell * P_sell(:,s)) + ... % 售电收益 sum(penalty * P_curtail(:,s))); % 弃风弃光惩罚 end % 约束 Constraints = []; for s = 1:ns for t = 2:nT % 功率平衡（并网微电网） Constraints = [Constraints, ... P_pv(t,s) + P_wind(t,s) + sum(P_g(t,:,s)) + P_bat(t,s) + P_ev(t,s) + P_buy(t,s) ... == P_load(t,s) + P_sell(t,s)]; % 储能SOC递推 Constraints = [Constraints, ... E_bat(t,s) == E_bat(t-1,s) + eta_bat * P_bat(t,s)]; % 燃气机组出力上下限 Constraints = [Constraints, ... P_g_min * x_unit(t,:) <= P_g(t,:,s) <= P_g_max * x_unit(t,:)]; % EV聚合功率上下限 Constraints = [Constraints, ... P_ev_min(t,s) <= P_ev(t,s) <= P_ev_max(t,s)]; end end % 求解 optimize(Constraints, objective, sdpsettings('solver','gurobi','verbose',2,'mipgap',0.01));

这里提醒几个细节。第一，YALMIP定义带场景下标的三维变量时内存涨得很快，场景数加到50以上、时段24、设备数若干，模型规模会到几千个变量，这时候尽量用矩阵化写法，避免在循环里反复用方括号拼接约束。第二，EV聚合功率P_ev我习惯写成“正为充电、负为放电”，约束里直接用P_ev_min ≤ P_ev ≤ P_ev_max限制，比分别定义充放电两个变量更简洁，但代价是目标函数里要加上EV电池损耗惩罚项，否则优化器会在低谷电价时让EV无节制充电、高峰时猛放电，场景循环往复，策略看似省钱实则毁电池。

4.3 求解器选型与常见线性化技巧

随机优化调度在大部分情况下会被建模成混合整数线性规划MILP或纯线性规划LP。涉及机组启停就必然有0/1变量，属于MILP。我一般用YALMIP做建模层，求解器优先选Gurobi或CPLEX；如果没有商业求解器，Matlab自带的intlinprog也能跑小规模案例，但超过一定规模会非常吃力。

用MILP求解就要注意线性化，最常见的是三类问题。第一，EV不能同时充电和放电。如果采用单变量P_ev的正负表示方向，这个问题天然规避；如果非要分成P_charge和P_discharge两个变量，就必须加互斥约束：

P_charge(t) ≤ M * z(t) P_discharge(t) ≤ M * (1 - z(t))

其中z是0/1变量，M取足够大的数，这叫大M法。M不能设得特别大，太大会造成数值病态，建议取该节点最大可能的功率交换量。

第二，购售电不能同时发生。并网微电网里要设定联络线最大交换功率，加入互斥约束：P_buy ≤ M * z_buy，P_sell ≤ M * (1 - z_buy)。这样避免模型出现“低价买入、高价卖出”同时存在的套路漏洞。

第三，二次型成本函数。燃气轮机燃料成本严格说应该是二次函数，但如果直接放进模型，MILP就变成MIQP，求解难度大幅上升。我的做法是分段线性化：把燃料成本曲线按出力区间切成3到5段，每段用斜率近似，这样模型保持MILP，求解稳定性和速度都好很多。Gurobi对MIQP的支撑虽然不错，但问题规模一大照样卡住，能线性化就线性化。

求解器参数也要会调。sdpsettings里设置mipgap=0.01，即容忍1%的次优性，往往就能把求解时间从几分钟压到几十秒。对工程调度来说，1%的误差完全可接受。