电动汽车接入区域综合能源系统之后原来教科书里那套“负荷曲线加机组组合做一个经济调度”的玩法就开始失灵了。充电桩扎堆、V2G反向送电、车主随机出行导致的负荷脉冲让配电网侧的净负荷曲线变得极难预测而区域综合能源系统又牵扯电、气、热、冷多条能量流的耦合优化。这篇博文是我把“含电动汽车的区域综合能源系统优化调度”从数学建模到Matlab实现完整走一遍的记录内容包括问题建模、求解器选型、Yalmip代码框架、算例结果分析以及我在实际编程里踩过的几个比较隐蔽的坑。适合正在做综合能源优化调度方向的研究生、做园区或配网级能量管理的工程师以及准备用Matlab做相关仿真的入门同学参考。1. 含电动汽车的区域综合能源系统问题本质与建模难点1.1 为什么电动汽车会让传统调度模型失效传统电力系统经济调度的核心假设是负荷可预测、电源可控顶多再加上风光的波动性。但电动汽车大规模接入后这个假设被撕开了两个口子。第一个口子是充电负荷的时空随机性。一辆EV什么时候充、在哪里充、充多久取决于车主的通勤习惯、剩余电量、充电桩空闲状态甚至当天的天气。单台车看着没多少容量但一个区域内几百上千台车的充电行为叠加起来就是一条锯齿状的负荷曲线早晚高峰期间的充电需求完全可以和原有的负荷高峰重叠导致配变过载。第二个口子是V2G带来的双向功率流动。电动汽车不再只是消费者它还可以在电价高峰时段向电网放电。这意味着调度模型里必须引入“充电”和“放电”两个方向的可控变量还要考虑电池循环损耗、车主参与意愿、离网时最低电量需求等约束。传统机组组合模型里的决策变量是连续的机组出力而EV充放电天然带有逻辑切换属性——同一时刻不能既充又放这就要引入二进制变量。如果把问题放大到区域综合能源系统情况更复杂电网、天然气网、热网通过燃气轮机、电锅炉、热电联产机组耦合在一起EV充电负荷的变化会通过电力平衡传导到燃气轮机的出力进而影响热功率输出和储热装置的动作。也就是说电动汽车扰动的不只是电力系统而是整个多能互补系统的运行状态。1.2 区域综合能源系统的典型拓扑结构与能量流分析我在搭建模型时采用了一个比较经典的结构也是文献里最常见的配置区域内有风电机组、光伏阵列、燃气轮机、电锅炉、蓄电池储能、蓄热罐、电动汽车聚合商以及不可控的电负荷和热负荷。外部电网通过联络线向区域供电天然气网络向燃气轮机供气。这个结构的能量流可以拆成三条主线电功率流风电、光伏、燃气轮机发电、储能放电、EV放电、电网购电共同为电负荷、电锅炉、储能充电、EV充电供电。热功率流燃气轮机余热回收、电锅炉制热共同为热负荷和蓄热罐供能。天然气流外购天然气供给燃气轮机转化为电和热两种二次能源。写到这里你应该能看出一个问题电动汽车的充放电策略会同时影响电力平衡和热力平衡因为它改变了燃气轮机的出力点。比如在夜间风电大发时段如果EV集中充电可以消纳一部分风电减少弃风而在傍晚负荷高峰时段如果EV放电可以减少从电网购电量但燃气轮机出力也会因此下降余热回收减少如果此时热负荷还在高位电锅炉就要补上热功率缺口。这种电-热耦合效应是单纯研究配电网EV调度时看不到的也是做区域综合能源系统调度最需要仔细处理的点。2. 优化调度数学模型目标函数与约束条件设计2.1 目标函数的选择单目标成本最小化是主流在做研究型仿真时最常见的目标是系统总运行成本最小化我这次也采用这个目标。之所以不推荐一开始就上多目标是因为EV调度涉及大量0-1变量和时序约束多目标优化比如成本碳排放双目标用Matlab实现时Pareto前沿的求解和决策者偏好设置会显著增加代码复杂度。先跑通单目标再扩展多目标是更稳妥的路线。我们的目标函数可以拆成四个部分购电成本从外部电网购买电量的费用一般用分时电价乘以联络线功率。购气成本燃气轮机消耗的天然气费用。设备运维成本风电、光伏、燃气轮机、储能、电锅炉的运维成本通常按功率或电量线性折算。弃风弃光惩罚为了保证新能源消纳在目标函数中加入弃风弃光量的惩罚项避免模型为了省钱而随意切掉风光出力。还有一个容易被忽略的部分是EV用户的参与补偿成本。如果调度中心要求EV在特定时段放电则需要支付给车主一定的激励费用否则车主没有动力参与V2G。这个成本可以直接加在目标函数里也可以在约束条件里设置强制放电量上限来间接体现。我的做法是把它按放电电量线性计入这样模型能够自行权衡系统是否值得为了一度电的峰谷差而付出激励成本。2.2 关键约束条件从电功率平衡到EV充放电边界约束条件是这个模型的重头戏每一个约束的物理意义都要搞清楚否则代码跑出来的结果可能就是“数学上可行、物理上荒谬”。电功率平衡约束是最基本的P_GT P_WT P_PV P_dis P_ev_dis P_grid P_load P_ch P_EB P_ev_ch这里P_GT是燃气轮机出力P_WT和P_PV是风电光伏实际出力P_dis和P_ch是储能放电和充电功率P_ev_dis和P_ev_ch是EV聚合体的放电和充电功率P_grid是联络线功率P_load是常规电负荷P_EB是电锅炉耗电。这个约束看起来简单但要注意单位统一我见过不少同学把kW和MW混在一起结果求解出来的功率差了三个数量级。储能SOC约束不能用简单的线性等式表达因为充放电效率和SOC上下限都要考虑SOC(t1) SOC(t) - P_ch(t) * eta_ch / E_cap P_dis(t) / (eta_dis * E_cap) SOC_min SOC(t) SOC_max还需要约束储能不能同时充放电这里要引入0-1变量u_ch和u_dis0 P_ch(t) u_ch(t) * P_ch_max 0 P_dis(t) u_dis(t) * P_dis_max u_ch(t) u_dis(t) 1电动汽车约束是本文的重点。EV不是单纯的储能它有出行需求。调度模型里必须设置每个EV的接入时段、初始SOC、目标SOC、最大充放电功率和电池容量。最关键的约束是EV在离开节点时SOC必须达到用户期望值。SOC_ev(k, t_dep) SOC_ev_target(k)这个约束如果漏掉模型会把EV的电池当成免费储能来用——白天电价高时就放电反正不担心车主开不走车。另外EV的到达时间、离开时间、初始SOC都来自出行链模拟需要提前生成好数据文件。燃气轮机和电锅炉约束包括出力上下限、爬坡约束P_GT_min P_GT(t) P_GT_max -P_ramp_GT P_GT(t1) - P_GT(t) P_ramp_GT 0 P_EB(t) P_EB_max有些文献还会加入最小启停时间约束但我在基础模型中先不做因为0-1变量一旦和爬坡约束耦合求解时间会大幅上升。先跑通基础模型再逐步加约束是仿真项目比较务实的推进方式。联络线功率约束限制了区域和大电网之间的交换功率-P_grid_max P_grid(t) P_grid_max这个约束体现了区域系统的自治性避免模型完全依赖大电网兜底。2.3 模型中的非线性项及其线性化处理上述模型看起来已经是线性的但实际写出来时会有几个隐藏的非线性点。最典型的例子是储能充放电效率的处理如果P_ch和P_dis是独立的连续变量那么效率可以直接乘在SOC递推公式里没有问题。但如果你用单一变量P_storage来表示储能功率正值放电、负值充电那么效率就得写成SOC(t1) SOC(t) - [P_storage(t) 0 ? P_storage(t) / eta_dis : P_storage(t) * eta_ch] / E_cap这个条件表达式没法直接放进线性规划框架。所以我在实现中采用P_ch和P_dis双变量加0-1标志位的方案虽然增加了变量数量但保证了模型的线性性求解稳定性高得多。另一个常被忽略的非线性来源是EV充放电功率和SOC的耦合。如果EV在充电过程中SOC达到上限模型必须允许停止充电而不是强制在接入时段内持续充电。因此EV充电功率的上限不能直接用固定值要写成分段形式0 P_ev_ch(k, t) u_ev_ch(k, t) * min(P_ev_ch_max, (SOC_ev_max - SOC_ev(k, t)) * E_cap_ev / eta_ch)这不是一个严格的线性约束但通过引入0-1变量并设置合理化上限可以用MILP方式求解。实际操作中我把EV的SOC上限约束单独写出再让充放电功率上限取固定值这样虽然略微保守但模型结构更干净。3. Matlab实现从数学建模到代码落地3.1 求解工具选型Yalmip加CPLEX还是启发式算法Matlab环境下做优化调度求解器的选型几乎决定了整个项目的走向。我的建议是优先使用Yalmip商用求解器CPLEX或Gurobi只在模型非线性程度极高时才考虑粒子群、遗传算法等启发式算法。原因有几点。第一Yalmip的建模语法非常接近数学表达写起来和推导公式几乎一一对应。第二CPLEX/Gurobi对MILP问题的求解速度和全局最优性保证是启发式算法无法比拟的。第三调试方便约束写错了可以打印出来逐条检查而粒子群跑完都不知道是模型错了还是参数没调好。如果你用的是MILP模型那直接上CPLEX。如果模型里包含二次项比如考虑网损的二次函数Gurobi处理二次约束更顺手。我这里用的是YalmipR2024aCPLEX 12.10的组合整体很稳定。注意CPLEX的许可证现在提供免费的学术版安装时选对版本否则在Matlab里调用时可能出现“指定路径无法找到Cplex”的提示。这个问题我遇到过不止一次基本都是因为只装了优化器本身没有把Matlab接口目录加到路径里。3.2 电动汽车充电负荷场景的生成方法EV优化调度的前提是有一组可信的EV行为数据。这部分我强烈建议用蒙特卡洛模拟生成而不是手动编造数据。具体步骤是统计区域内的EV数量、电池容量分布、充电功率等级。基于私家车出行规律假设EV每天有两个主要停放时段白天工作时段和夜间居家时段。对每台EV随机抽样到达时间、离开时间、初始SOC。我把这些数据保存成一个矩阵每一行代表一台EV列依次是电池容量、接入时段、离开时段、初始SOC、目标SOC、最大充放电功率。这样在Yalmip建模时可以直接按行索引变量代码非常清爽。%% 生成EV出行场景 % 输入EV数量、电池容量范围、通勤规律参数 N_ev 200; E_cap_ev 40 20 * rand(N_ev, 1); % 电池容量 40~60 kWh arrival_h 8 1.5 * randn(N_ev, 1); % 到达时段均值8点 depart_h 17 1.2 * randn(N_ev, 1); % 离开时段均值17点 soc_init 0.3 0.3 * rand(N_ev, 1); % 初始SOC 30%~60% soc_target 0.8 0.1 * rand(N_ev, 1); % 目标SOC 80%~90%这里生成的数据是连续值实际使用时需要把到达和离开时段折算成调度间隔的整数步长比如调度间隔为1小时那么8.3小时就按8处理。3.3 核心代码框架变量定义、约束构建与求解下面给出一段精简但可运行的Yalmip框架代码完整项目代码太长这里只展示建模核心逻辑。%% 定义优化变量 T 24; % 调度周期 % 常规设备 P_GT sdpvar(T, 1, full); % 燃气轮机出力 P_EB sdpvar(T, 1, full); % 电锅炉耗电 P_grid sdpvar(T, 1, full); % 联络线功率正为购电 % 储能 P_ch sdpvar(T, 1, full); % 充电功率 P_dis sdpvar(T, 1, full); % 放电功率 u_ch binvar(T, 1); % 充电标志 u_dis binvar(T, 1); % 放电标志 SOC sdpvar(T1, 1, full); % 储电SOC % 电动汽车每台EV接入时段内充放电功率 P_ev_ch sdpvar(T, N_ev, full); P_ev_dis sdpvar(T, N_ev, full); u_ev_ch binvar(T, N_ev, full); u_ev_dis binvar(T, N_ev, full); %% 目标函数 Objective sum(price_grid .* P_grid) ... sum(c_gas * P_GT) ... sum(c_om .* (P_GT P_WT P_PV)) ... penalty_curtail * sum(P_wt_max P_pv_max - P_WT - P_PV) ... c_ev_comp * sum(P_ev_dis(:)); %% 约束 Constraints []; % 电功率平衡 Constraints [Constraints, P_GT P_WT P_PV P_dis sum(P_ev_dis,2) P_grid ... P_load P_ch P_EB sum(P_ev_ch,2)]; % 储能约束 Constraints [Constraints, SOC(2:T1) SOC(1:T) - P_ch.*eta_ch/E_cap P_dis/(eta_dis*E_cap)]; Constraints [Constraints, SOC(1) SOC_init, SOC_min SOC SOC_max]; Constraints [Constraints, 0 P_ch u_ch*P_ch_max, 0 P_dis u_dis*P_dis_max]; Constraints [Constraints, u_ch u_dis 1]; % EV约束 for k 1:N_ev t_arr round(arrival_h(k) * 60 / dt) 1; t_dep round(depart_h(k) * 60 / dt) 1; % 仅接入时段内可充放电 for t 1:T if t t_arr || t t_dep Constraints [Constraints, P_ev_ch(t,k) 0, P_ev_dis(t,k) 0]; else Constraints [Constraints, 0 P_ev_ch(t,k) u_ev_ch(t,k)*P_ev_ch_max]; Constraints [Constraints, 0 P_ev_dis(t,k) u_ev_dis(t,k)*P_ev_dis_max]; Constraints [Constraints, u_ev_ch(t,k) u_ev_dis(t,k) 1]; end end % 离网时SOC达到目标值 SOC_ev_t soc_init(k) * E_cap_ev(k) ... sum(P_ev_ch(1:T,k).*eta_ch_ev - P_ev_dis(1:T,k)/eta_dis_ev) * dt / 60; Constraints [Constraints, SOC_ev_t soc_target(k) * E_cap_ev(k)]; end %% 求解 ops sdpsettings(solver, cplex, verbose, 2); result optimize(Constraints, Objective, ops);这段代码展示了整个求解框架。有几个细节需要特别说明dt如果按分钟调度公式里的转换关系就很重要我在EV离网SOC计算中乘了dt/60避免kWh和kW之间的单位错误。每台EV的接入时段约束写成了双重循环当EV数量很大时比如上千台约束数量会暴增求解时间可能从几秒涨到几分钟。这时候可以考虑把EV聚合成几个集群每个集群用聚合的功率和SOC描述代价是损失一些单台EV的SOC精度。3.4 结果输出与绘图求解完成后我通常会把结果整理成几个结构体方便后续统计和画图%% 结果整理 Result.P_GT value(P_GT); Result.P_grid value(P_grid); Result.P_ch value(P_ch); Result.P_dis value(P_dis); Result.SOC value(SOC); Result.P_ev_ch value(P_ev_ch); Result.P_ev_dis value(P_ev_dis); Result.Objective value(Objective);画图时我习惯用stairs画阶梯状的功率曲线因为调度结果在每个时段是恒定的用连续曲线反而会误导读者以为功率之间是平滑过渡的。SOC轨迹用plot加数据点标注就可以了。4. 算例测试与结果分析一个典型区域场景4.1 算例参数设置为了验证模型我用了一个典型区域的算例数据。区域内共有200台电动汽车风电装机2000kW光伏装机1000kW燃气轮机额定出力1500kW电锅炉额定功率800kW蓄电池容量1000kWh最大充放电功率200kW蓄热罐容量500kWh。外部电网联络线功率上限3000kW分时电价采用峰谷电价结构峰时8:00-11:00和18:00-21:00为1.2元/kWh平时6:00-8:00、11:00-18:00和21:00-23:00为0.7元/kWh谷时23:00-6:00为0.35元/kWh。燃气轮机效率取0.45热电比1.2天然气价格按每立方米折合到单位热值后计算。EV的电池容量40-60kWh最大充放电功率7kW充电效率0.95放电效率0.9车主目标SOC设为80%-90%。4.2 三种运行模式的对比结果我设置了三个对比场景场景A无序充电。EV接入后立刻以最大功率充电直到达到目标SOC不做任何优化调度。场景B有序充电。EV充电功率由优化模型调度但不允许放电。场景CV2G双向调度。EV既可以充电也可以放电参与系统削峰填谷。优化结果如下表所示场景总运行成本/元弃风弃光电量/kWh联络线高峰功率/kWEV放电量/kWhA无序充电52863116029500B有序充电4932862025800C V2G调度461353102210860从数值上能看出几个非常典型的结论第一有序充电相比无序充电可以降低约3500元的运行成本主要原因是把充电负荷转移到了谷时段减少了高峰期的购电费用。第二加入V2G放电后成本进一步降低约3200元放电量为860kWh这部分电量全部在18:00-21:00的高峰时段释放替代了部分高峰购电。第三弃风弃光率发生了显著变化无序充电场景下夜间谷时段EV充电功率有限风电大发时只能弃风而有序充电场景中模型会主动把EV充电安排到夜间相当于给风电提供了额外消纳空间。4.3 电-热耦合效应在算例中的表现我在数据分析时特别关注了一个细节V2G放电时燃气轮机的出力是否下降、热负荷缺口由谁填补。结果发现在傍晚EV集中放电的时段燃气轮机出力确实下调了约300kW余热回收减少后热负荷缺口由电锅炉增加了约220kW的耗电来弥补。也就是说EV放电节省的购电成本有一部分被电锅炉多耗的电抵消了。这就是区域综合能源系统和纯电力系统优化调度最本质的区别不能只看电力平衡要多能互补地看全局成本。如果我没有建立电-热耦合模型而只是简单地把EV当作配电网的灵活性资源就会高估V2G的经济效益。4.4 关键参数的敏感性分析我还做了一组敏感性分析重点观察EV数量对调度结果的影响。EV数量从100台逐步增加到500台每100台测算一次。可以看到随着EV数量增加总运行成本先下降后上升转折点大概在300台左右。原因是前期EV数量较少时EV作为灵活性资源的价值大于它的充电成本所以系统总成本下降但当EV数量超过一定阈值后额外增加的充电负荷对电网和燃气轮机形成的压力超过了灵活性收益成本反而上升。这个现象在文献里不太常被讨论但在实际项目里非常重要它直接影响一个区域该配置多少充电桩、是否值得推广V2G。提示如果你在做敏感性分析时发现EV越多越省钱大概率是模型里漏了EV充电的电能量成本或者漏了配电网容量约束。EV毕竟不是免费的储能它的每一次充电都要消耗能量只是转移了时间不会凭空增加系统效率。5. 我在实际Matlab编程中踩过的坑与调试经验5.1 Yalmip和求解器的版本匹配坑做这个项目时我踩过的第一个坑是Yalmip版本和CPLEX求解器的接口问题。一开始我装的是Yalmip 2022版本CPLEX是12.10在Matlab 2021b环境下跑结果是optimize返回一条“CPLEX interface not found”的错误。排查了半天发现新版Yalmip对CPLEX 12.10的支持已经标记为deprecated需要切换到Gurobi或者升级CPLEX。解决办法是把Yalmip更新到最新版或者在sdpsettings里显式指定求解器路径ops sdpsettings(solver, cplex, cplex.path, C:\Program Files\IBM\ILOG\CPLEX_Studio1210\cplex\bin\x64_win64);如果项目不强制要求CPLEX我更推荐直接使用Gurobi它和Yalmip的兼容性最好求解MILP的速度也比CPLEX快不少授权申请也比较方便。5.2 变量维度不一致导致的“隐形错误”Yalmip有个特点是矩阵运算自动识别维度但这也带来了一个隐蔽的问题如果你的某个变量维度是T而另一个是T1在构成约束时Yalmip不会报错而是自动广播成T1维导致最终约束数量比预期多了一行结果却依然能求解出来。最开始我在写SOC递推约束时就把SOC定义成了T1P_ch定义成T结果打印约束时发现SOC的初值约束没有生效因为被广播后的约束错位了。我的建议是在写完约束后用size(Constraints)验证维度或者用check(Constraints)逐条检查可行性。尤其当模型跑出来的最优成本明显不合理时比如为负第一反应应该是检查有没有隐式广播或漏了约束而不是怀疑求解器。5.3 大M法中的参数选取经验在建模储能充放电互斥约束时传统写法是u_ch u_dis 1但如果模型里还包含其他逻辑约束比如“储能SOC低于某值时禁止放电”就需要引入大M法P_dis(t) M * u_dis(t)这里的M如果取得太大比如取1e6会严重恶化MILP的松弛性能导致求解时间从几秒暴涨到几分钟如果取得太小又可能把可行域切掉一小块得到局部最优。我的经验是M取该变量物理意义上限的2到3倍即可比如EV放电功率上限是7kWM取20就够了不要为了保险而盲目放大。仿真中凡是遇到M参数都应该做一次敏感性测试看看结果是否对M的取值敏感。5.4 初始SOC设定对EV调度结果的干扰EV模型的初始SOC如果设得太高比如适应度测试时随手设成90%那么EV在调度周期内几乎不需要充电V2G放电潜力也极其有限得到的结果就是“EV对系统无影响”这显然是错的。反过来如果初始SOC设得太低比如10%模型会强迫EV尽快充电可能在前几个时段产生一个充电脉冲破坏削峰填谷的效果。合理的做法是让初始SOC服从一个分布比如我前面提到的30%-60%均匀分布并让目标SOC在80%-90%之间随机分布。这样单台EV的行为有差异性聚合起来的充电需求才真实。5.5 求解时间优化从40分钟到2分钟的改造这个项目最初跑一个500台EV的场景需要超过40分钟后来优化到2分钟以内。主要的优化手段有三个用EV聚合模型替代逐台建模将200台EV按照到达时段和离开时段聚类成5-8个集群每个集群内用平均SOC和总充放电功率表示将变量数量减少一个数量级。去掉非必要的0-1变量。比如储能互斥约束如果不影响实际最优解可以尝试用互补旋转约束替代或者直接删掉互斥约束靠目标函数的价格差自然避免同时充放电。调整求解器的MIP gap。学术研究中强求gap0%确实能拿到严格最优但实际项目里gap0.5%和0在成本上的差异几乎可以忽略不计。我把Gurobi的MIPGap设为0.005后求解速度提升了将近8倍。我在实际项目中体会很深的一点是做这类优化调度研究模型跑通只是第一步后面大量的时间会花在调试约束、检查数据、验证结果合理性上。而严谨的调试习惯比如每次修改模型后先跑一个小规模算例再上完整规模能替自己省下一整天的排查时间。这个“先小后大”的策略是学长教我的现在轮到我把它在这里分享出来。另外如果你准备沿着这个方向继续扩展比较推荐的下一步是加入碳交易机制或者多目标优化成本与碳排放甚至把EV的出行链和交通路网耦合进来那个问题就更有意思了。