一个不算太复杂的标题背后其实藏着一整套需要捋顺的逻辑链条新能源出力随机性怎么进模型、电热气多能流怎么耦合、优化目标怎么定、最后怎么用 Matlab 把这套东西算出来。我当初拿到这个题目的时候第一反应是又是个披着综合能源外壳的迭代求解器但真正下手拆完才发现问题的核心不在算法多花哨而在不确定性处理方式的选择和多能流耦合约束的建模精度上。这篇文章我就按我实际做项目的顺序来写把我踩过的坑、验证过的做法都交代清楚想复现这个项目的朋友可以直接按这个思路走。1. 项目定位为什么综合能源系统优化不能按电网调度的老思路做综合能源系统Integrated Energy SystemIES和传统电力系统调度的本质区别在于能量不只有电一路。一个典型的园区级综合能源系统里电网、天然气网、热网是同时存在并且通过耦合设备绑在一起的。电锅炉、热泵、燃气轮机、热电联产机组CHP、电储能、蓄热罐这些电气设备就是不同能源网络之间的连接器和转换器。如果只优化电网侧那就漏掉了热负荷对电负荷的反向约束——比如冬天采暖需求上来了电锅炉一开电负荷曲线立刻被拉出一个高峰。反过来如果只优化热网侧燃气锅炉的耗气量又会影响整个系统的运行成本。所以协同优化要做的是把这个三角关系放到一个统一框架里同时决策每个时刻每台设备出多少力、储能充放多少、从电网买多少电、从气网购多少气。这个项目标题里还有两个关键词值得注意不确定性和协同。不确定性是指风电和光伏的出力不是定值预测曲线和实际曲线之间存在偏差如果优化模型把预测值当确定量处理调度计划在实时运行时会存在功率失衡风险。协同则是指电网、气网、热网之间的耦合约束要同时满足而不是各自独立优化后再去凑结果。所以这个项目的实际目标可以拆成三块建立一个含多种电气设备的多能流耦合模型覆盖电、热、气三种能量形式把新能源出力不确定性以数学上可处理的形式嵌入优化模型以系统总运行成本最小为目标求出一个满足所有设备出力约束和网络平衡约束的日前调度方案并输出各设备的逐时出力曲线。适合参考这个项目的读者我大致分三类一是做综合能源系统调度方向的研究生需要一套能跑通、能出图的基线代码二是做园区能源管理的工程人员想了解源网荷储一体化的优化模型怎么落地三是准备用 Matlab 做能源系统优化但还没摸清 Yalmip 建模技巧的初学者。三类读者的切入点不同但最终都要回归到同一个核心问题上模型里的约束和变量怎么组织才能让求解器高效地算出来。2. 新能源出力不确定性建模从点预测到随机场景2.1 确定性模型的问题在哪如果不考虑不确定性风电场和光伏电站的出力在模型里就是一个固定的时间序列比如24小时的逐时预测功率曲线。这种做法最简单优化出来的结果也看起来最优,因为它在数学上只面对一个确定性的场景。但在实际运行中光伏出力受云层遮挡影响可能某天上午的真实出力比预测值低30%风电更是如此风速预测偏差大出力波动范围经常超过50%。一旦真实出力低于计划值系统就面临两类后果一是供给侧的缺口要靠快速备用补上而备用成本往往远高于正常调度成本二是如果备用不足只能切除负荷这在工程上是不可接受的。所以不考虑不确定性的调度方案本质上是拿一个理想剧本去指导一个有随机扰动参与的实际过程两者之间的差距就是弃风弃光风险和失负荷风险。不确定性建模要解决的核心问题就是把这个差距提前码进优化模型里让调度方案自带鲁棒性。2.2 常见的三类不确定性建模方法处理新能源出力不确定性业界和学术界常用的方法主要有三类我在做这个项目时也对比过它们的适用场景。场景法Scenario-based。用一组带有概率权重的出力场景来代替单条预测曲线。场景来源于历史出力数据的统计聚类或者是对预测误差分布进行采样。优点是直观模型规模可控求解难度适中缺点是场景数量直接影响计算量——场景太少保不住多样性场景太多求解时间线性上涨。鲁棒优化Robust Optimization。不关心具体概率分布只关心最坏情况。给定一个不确定性区间比如光伏出力波动不超过预测值的±20%模型在这个区间内寻找最坏情况下的最优解。优点是非常保守可靠缺点也正是过于保守——按极端场景做的调度计划在大多数正常日子里的运行成本会偏高。区间预测 机会约束Interval Forecast Chance Constraint。通过预测模型给出每个时刻出力的置信区间再把区间不被破坏的要求转化为概率约束允许以小概率违反约束来换取更低的经济成本。这条路径处理起来相对复杂对求解器的支持要求也更高。我在这个项目里采用的是场景法加缩减处理的组合套路。原因很直接在 Matlab 环境下场景法可以与 Yalmip 建模框架无缝对接配合 CPLEX 或 Gurobi 求解器可以稳定拿到全局最优解代码也好调试。而且场景缩减技术比如基于Kantorovich距离的同步回代缩减可以先把1000个场景缩减到20个左右既保留出力波动的多样性又控制模型规模。2.3 场景生成与缩减的工程实现场景生成的做法是基于预测误差的概率分布进行蒙特卡洛采样。假设某个时刻风电预测出力是 (P_{w}^{pre}(t))预测误差 (\varepsilon(t)) 服从均值为零的正态分布标准差取该时刻预测出力的一个比例比如15%那么采样场景 (s) 在 (t) 时刻的风电出力就是[ P_{w,s}(t) \max\left(0,; P_{w}^{pre}(t) \varepsilon_s(t)\right) ]这里要特别强调一下采样必须逐时刻进行且要加入时间相关性。如果每个时刻独立采样生成出的场景在时间轴上会出现锯齿状跳变整个场景的出力曲线看起来像白噪声既不符合实际风速变化的连续性也会让优化结果变得不可理解。处理方式是在误差序列上做一阶自回归处理让相邻时刻的误差保持相关性[ \varepsilon_s(t) \rho \cdot \varepsilon_s(t-1) \eta_s(t) ]其中 (\rho) 是相关系数我一般取0.8到0.9(\eta_s(t)) 是新的随机扰动项。这样生成出来的场景就有了一天之内平滑波动的形态更接近真实的出力过程。场景缩减环节最简单有效的是快速前向选择法。基本逻辑是从初始场景集合里每一步选中一个与其他场景概率距离之和最小的场景作为代表场景然后更新剩余场景的概率权重。我用 Matlab 写这个缩减算法时一般控制目标场景数为初始场景数的5%左右比如800个场景缩到40个。缩减后的场景不仅能覆盖大部分极端情况而且概率分布和大规模场景集合的距离也控制在可接受范围内。这个预处理环节虽然不直接产生优化结果但它决定了优化结果对不确定性描述得够不够真实值得多花点时间做扎实。3. 协同优化模型搭建目标函数与约束的关键取舍3.1 优化目标总运行成本到底包含哪些项模型的目标函数我在项目里定义的是系统日运行总成本最小化它的构成不是单纯几台设备的运行费用相加而是要覆盖以下几个层面购电成本从上级电网购入电力的费用按分时电价计算这是成本的绝对大头购气成本从气网购入天然气的费用按单位热值价格折算设备运行维护成本各种电气设备在出力过程中产生的可变运维费通常按出力大小的比例系数计算弃风弃光惩罚成本目的在于防止优化结果出现宁买电也不用免费新能源的反常行为失负荷惩罚成本给系统加了保险丝防止模型为了省钱而把负荷削掉。目标函数用数学形式写出来是[ \min ; \sum_{t1}^{T} \left( c_e(t)P_{grid}(t) c_g F_{gas}(t) \sum_{i} c_{om,i} P_i(t) c_{cut} P_{curtail}(t) c_{loss} P_{loss}(t) \right) ]我在这里重点说一个实践中容易被忽略的细节设备运维成本系数不能拍脑袋定。很多人建模时随便给个0.01元/kWh结果优化结果里设备频繁启停因为启停不产生成本模型自然不会约束。真实场景中设备启停会带来额外的磨损和维护费用所以我在目标函数里通常会加一项启停惩罚项让模型自动避开无意义的频繁切换。3.2 核心约束多能流耦合约束是协同优化的灵魂约束条件按作用范围我习惯把它拆成四组。每一组约束的背后都对应着一个物理过程或者运行边界。第一组是电力平衡约束。每个时刻系统的电功率必须守恒。等式约束写成[ P_{grid}(t) P_{wind,s}(t) P_{pv,s}(t) P_{chp,e}(t) P_{dis}(t) L_e(t) P_{ch}(t) P_{eb}(t) P_{hp}(t) P_{curtail}(t) ]从左到右依次是电网购电、风电、光伏、CHP电出力、储能放电右侧是电负荷、储能充电、电锅炉耗电、热泵耗电和弃电。这里特别注意对于每个新能源出力场景 (s)这个等式都要满足所以这个约束的真实数量是时刻数乘以场景数。第二组是热力平衡约束。热负荷由CHP的余热、燃气锅炉、电锅炉、热泵共同满足蓄热罐可以充当缓冲[ H_{chp}(t) H_{gb}(t) H_{eb}(t) H_{hp}(t) H_{tank,dis}(t) L_h(t) H_{tank,ch}(t) ]热力平衡约束的难点在于热负荷本身有延迟性如果模型里没有储热设备那么热电联产机组的以热定电特性会直接把电出力也锁死导致优化空间被压缩。加蓄热罐的意义就是解耦热电的强耦合关系让系统在电价低谷时产热储热、电价高峰时放热保供电。第三组是设备出力约束。每一台设备都有上下限约束[ P_i^{min} \le P_i(t) \le P_i^{max} ]以及爬坡约束[ -\Delta P_i^{down} \le P_i(t) - P_i(t-1) \le \Delta P_i^{up} ]风电和光伏的出力上限还要进一步受制于该场景下的可用出力当新能源出力多了、系统无法完全消纳时才启用弃风弃光变量。第四组是储能约束包括储能SOC状态转移方程[ E_{store}(t1) E_{store}(t) \eta_{ch} P_{ch}(t) - \frac{P_{dis}(t)}{\eta_{dis}} ]以及容量上下限约束、充放电功率上限约束。这里有个建模上的坑就是充放电不能同时进行——如果不加约束模型会因为储能效率都小于1而天然不选择同时充放可以用一组二进制变量来严格约束也可以直接用效率逻辑把这个行为写死。3.3 与确定性模型的对比方案我在项目里额外做了一套对比实验一套是考虑不确定性的场景随机优化一套是把新能源出力固定为预测值的确定性优化。这样对比的好处是可以直接量化不确定性对调度成本和调度方案的影响。从结果上看确定性模型算出来的日运行成本确实更低但它在真实出力发生偏差时系统需要额外调用备用资源实际的运行成本远高于名义成本。而场景随机模型因为提前考虑了多种可能性算出来的调度方案虽然名义成本高一些但方案在真实场景下执行时的后悔值更小。如果系统本身有储能和蓄热这类灵活资源随机模型的这一优势会被放大因为优化过程会自动把多余新能源出力导向储能设备为可能发生的低谷时段提前蓄力。4. Matlab代码实现从数学公式到可运行程序4.1 程序框架和核心模块整个项目的 Matlab 代码我是按照模块化思路组织的文件结构大致如下main.m主程序设置基础参数、调用各模块、汇总结果并绘图load_data.m加载负荷曲线、新能源预测出力、分时电价和天然气价格等数据文件scenario_generation.m蒙特卡洛采样生成新能源出力场景并执行场景缩减build_model.m构建 Yalmip 优化模型定义变量、目标函数和约束条件solve_model.m调用 CPLEX 求解器求解并处理求解状态plot_results.m输出各设备出力曲线、电网交互功率、储能SOC和系统总成本等图表。主程序里跟日期相关的参数全部集中配置在开头包括优化周期、时间分辨率我用的1小时、负荷数据文件的路径和场景数量等。这样后期如果要改成15分钟分辨率只需要改时间参数模型会自动扩展约束数量不需要改动其他逻辑。4.2 Yalmip建模的关键写法Yalmip 是 Matlab 环境下一个非常好用的建模工具底层对接 CPLEX、Gurobi 等多个求解器。它的价值在于把优化变量的定义、约束的组装和求解器的调用统一封装起来代码可读性比直接调 CPLEX API 高很多。定义变量的时候要注意区分连续变量和二进制变量。设备出力、储能SOC这些是连续变量储能充放电状态、设备启停状态是二进制变量。具体到项目里我这样写% 定义优化变量 % P_chp: N_scenarios x T 的矩阵表示各场景下CHP电出力 % P_gb: 燃气锅炉热出力 % P_eb: 电锅炉耗电功率 % P_hp: 热泵耗电功率 % P_w_use: 各场景下实际消纳的风电功率 % P_pv_use: 各场景下实际消纳的光伏功率 % P_grid: 从电网购电功率 % u_chp: 二进制变量CHP启停状态 P_chp sdpvar(N_scenarios, T, full); P_gb sdpvar(1, T, full); P_eb sdpvar(1, T, full); P_hp sdpvar(1, T, full); P_w_use sdpvar(N_scenarios, T, full); P_pv_use sdpvar(N_scenarios, T, full); P_grid sdpvar(1, T, full); u_chp binvar(1, T, full);这里的一个关键点是凡是受新能源场景影响的变量都定义为场景数x时段数的矩阵不受场景影响的决策变量电网购电、设备出力则定义为1xT的向量。这种做法在阶段优化的框架里就是所谓的第一阶段的决策变量对所有场景统一第二阶段的变量按场景分别取不同值。如果把这个结构搞混了模型求解结果会出现一个场景下的决策偷看了另一个场景的信息这本质上是非预期性问题会让优化结果失真。约束条件的写入用 Yalmip 的循环加赋值方式。以一个关键约束为例电力平衡约束可以写成Constraints []; for s 1:N_scenarios Constraints [Constraints, ... P_grid P_w_use(s,:) P_pv_use(s,:) P_chp(s,:) * eta_chp_e P_dis ... L_e P_ch P_eb P_hp P_curtail(s,:)]; end特别注意像电网购电功率 (P_grid) 这种不依赖场景的变量在循环里直接索引出来的效果是和每个场景都关联这正好满足约束要求而如果错误地把 (P_grid) 也定义成场景矩阵就会得到每个场景不同的购电决策那就偏离了协同优化的初衷。4.3 求解配置与参数调试模型组装好之后调用求解器的代码非常简洁options sdpsettings(solver, cplex, verbose, 2, showprogress, 1); optimize(Constraints, Objective, options);这里常见的一个问题是求解时间不可控。全场景不缩减、直接塞给 CPLEX 的时候模型规模可能达到几万约束和几千个变量其中还包含二进制变量求解时间可能从几十秒暴涨到几十分钟。我实测的经验是场景数量控制在20到30个之间二进制变量只用于储热罐和启停状态数量级几十个CPLEX 通常能在几十秒内完成求解。如果遇到求解时间过长优先检查的不是求解器而是模型里二进制变量的数量是不是膨胀了。比如储能不能同时充放这个约束可以用两个二进制变量的相互约束实现但如果每个场景每时刻都引入的话二进制变量规模就直接翻倍求解时间会呈指数级上涨。我处理这类约束的做法是由于所有场景共享一组储能状态二进制变量只在一个代表场景里定义充电状态用逻辑约束关联到其他场景这样能有效控制整数变量的规模。4.4 数据的组织方式这个项目对数据格式的要求是必须提前明确好的。负荷数据、风电预测数据、光伏预测数据、分时电价和天然气价格要统一组织成列向量或矩阵时区必须对齐。我踩过一次亏负荷数据的第一个点是0点而新能源预测数据从第1小时开始记录两者错位了1个小时导致优化结果第一小时和最后一小时出现完全不符合物理的跳变。排查了半天才发现是数据索引错位。数据文件我倾向于用.mat格式存储在load_data.m里一次性加载并立即做单位统一。比如电价用元/kWh天然气价格用元/m³但CHP的燃气耗量是m³热值转换是每立方米天然气的热值对应多少kWh这些单位如果不提前统一目标函数里的成本项相加就会彻底乱掉。5. 复现跑通后的验证思路与实操踩坑记录5.1 怎么确认结果是对的拿到优化结果之后第一件事不是看成本降低了多少而是做物理合理性检验。我总结了四个必须检查的项能量平衡验证把优化结果中的各项代回电力平衡和热力平衡方程检查残差是否为零。Yalmip 可以直接对变量赋值后计算约束残差这一步能在五分钟内验证约束是否真的被满足。出力边界检查所有设备的出力是否符合上下限约束。如果某台设备的出力曲线平滑地压在最大值线上这是正常的如果出现超出边界的点那一定是建模问题。储能SOC连续性检查蓄热罐的SOC曲线是否出现阶梯性跳变。理论上SOC变化应当符合能量转移方程跳变意味着SOC约束写得不对。弃风弃光行为分析在新能源出力特别大的时段模型是否优先消纳新能源。如果弃风弃光惩罚成本设置合理只有在极端情况下才会出现弃电行为。5.2 几个典型坑我复现这个项目的时候遇到过几个坑写出来给大家省时间。第一个坑是蓄热罐的状态变量初始值没给。Yalmip 建模的时候储能SOC的递推约束第一时刻需要初始SOC参数漏了这个参数会导致第一时刻的约束残缺求解结果跟后续时刻的SOC断层。解决方式是提前定义E_store_init参数并把它作为常数项放进约束等式右边。第二个坑是目标函数里的成本项量纲不一致。CHP的燃气消耗量如果用m³表示而电价用元/kWh 表示那么最终目标函数就是元 元看似没有报错实际是把燃气成本单位搞混了。应该统一换算成 kWh 热值或者统一按能量单位计费再去相加。第三个坑是场景缩减后概率权重没有归一化。缩减算法的目的是从原始场景集中选出代表性场景并重新分配概率如果代码里少了归一化步骤目标函数会因为没有概率权重或权重和为1导致成本计算结果整体偏小或偏大和真实期望成本完全对不上。第四个坑是优化结果里出现新能源出力和用电负荷同时出现尖峰。排查后发现问题出在负荷预测数据里包含了某些异常值优化模型为了满足尖峰负荷被迫让所有机组满出力结果成本和出力曲线都变得不合理。最终是做了数据清洗把明显的异常点做平滑处理之后结果才正常。这里也建议大家拿到负荷数据的第一件事就是画一遍时序曲线用眼睛先做一轮质检。5.3 结果对比的正确打开方式项目里我对比了三组实验结果确定性优化、场景随机优化不考虑储能、场景随机优化考虑储能。三组结果并排对比从数据上能看到三条核心规律考虑不确定性之后电网购电曲线的峰谷形态会发生变化模型更倾向于在电价低时段多买电并储存起来用来对冲新能源出力的随机波动储能设备的引入可以把弃风弃光率下降5到10个百分点同时系统总成本降低约3%到5%具体数值取决于储能容量和目标函数的惩罚项设置场景数从10个增加到30个时成本结果逐渐收敛从30个继续增加结果变化很小但计算时间增长明显。这也可以作为大家调参时的经验参考——场景数没必要盲目往大里加够用就行。6. 进一步扩展的方向代码跑通之后可以走的方向就比较开阔了。以我自己的实践来看最有价值的扩展有三个。一是引入碳交易成本。当前目标函数只考虑了经济成本如果把碳排放配额、碳交易价格纳进去系统的调度决策会发生明显变化。天然气单价高但碳排放相对低的机组在碳约束下会更有竞争力而煤电占比高的购电来源会变得不再经济整个调度方案的成本结构会被重新洗牌。二是将时间分辨率从1小时提升到15分钟。更细的时间分辨率可以让储能设备的充放电策略更精确减少跨时段功率缺口问题。模型约束数会翻四倍但现代求解器和Yalmip框架处理这种规模没有太大压力。我建议在代码把时间分辨率做成参数一劳永逸。三是引入多园区之间协同。现在的模型是单园区独立优化如果扩展为多园区的综合能源系统协同优化就需要加入网络拓扑约束和共享储能设备模型约束的数量会显著增加但实际的工程价值会更大因为多园区之间的能源互补潜力通常比单园区内部的优化空间更可观。就我自己实测下来的体会这类项目的核心难点从来不在求解器调用和代码语法上而是在模型假设的合理性上。新能源出力的不确定性、设备效率的参数设置、负荷数据的质量——每一个环节的细节偏差都会在层层叠加之后放大为整体结果的失真。写代码之前一定要先把物理过程在脑子里模拟一遍让每一个变量和约束都落到实处后面跑通了才能不退回去返工。